AD-A261li608
of Engineers
DREDGING RESEARCH PROGRAM
TECHNICAL REPORT DRP-92-6
ADCIRC: AN ADVANCED THREE-DIMENSIONAL
CIRCULATION MODEL FOR SHELVES,
COASTS, AND ESTUARIES
Report 1
THEORY AND METHODOLOGY
TIC
-ECTE I
.R3 19931
OF ADCiRC-2DDI AND ADCIRC-3DL
R. A. Luettich, Jr.
Institute of Marine Sciences
University of North Carolina at Chapel Hill
Morehead City, North Carolina 27514
J. J. Westerink
Department of Civil Engineering
University of Notre Dame, Notre Dame, Indiana 46556
Norman W. Scheffner
Coastal Engineering Research Center
DEPARTMENT OF THE ARMY
Waterways Experiment Station, Corps of Engineers
3909 Halls Ferry Road, Vicksburg, Mississippi 39180-6199
93-04407
iiiliil iin
November 1992
Report 1 of a Series
Approved For Public Release; Distribution Is Unlimited
Prepared for DEPARTMENT OF THE ARMY
US Army Corps of Engineers
Washington, DC 20314-1000
Under Work Unit No. 32466
98 82 078
The Dredging Research Program (DRP) is a seven-year program of the US Army Corps of Engineers.
DRP research is managed in these five technicat areas:
Area 1 - Analysis of Dredged Material Placed in Open Water
Area 2 - Material Properties Related to Navigation and Dredging
Area 3 - Dredge Plant Equipment and Systems Processes
Area 4 - Vessel Positioning. Survey Controls, and Dredge Monitoring Systems
Area 5 - Management of Dredging Projects
Destroy this report when no longer needed. Do not return
it to the originator.
The contents of this report are not to be used for
advertising, publication, or promotional purposes.
Citation of trade names does not constitute an official
endorsement or approval of the use of such
commercial products.
Dredging Research Program
Report Summary
of Engineers
Waterways Experiment
Station
US Army Corps
ADCIRC: An Advanced Three-Dimensional Circulation Model for Shelves, Coasts,
and Estuaries; Report I, Theory and Methodology of ADCIRC-2DDI and
ADCIRC-3DL (TR DRP-92-6)
ISSUE: A unified and systematic methodol¬
ogy must be provided to use in the investiga¬
tion of the dispersive or nondispersive charac¬
teristics of a site proposed for the disposal of
dredged material in open water as well as to
analyze existing disposal sites.
RESEARCH: ADCIRC (Advanced Three-
Dimensional Circulation Model) was devel¬
oped as a part of the Dredging Research Pro¬
gram (DRP) as a means of generating a
database of harmonic constituents for tidal ele¬
vation and current at discrete locations along
the east, west, and Gulf of Mexico coasts and
to utilize tropical and extratropical global
boundary conditions to compute frequency-
indexed storm surge hydrographs along the
US coasts. The database is being developed
to provide site-specific hydrodynamic bound¬
ary conditions for use in analyzing the long¬
term stability of existing or proposed dredged
material disposal sites.
SUMMARY: The report describes the the¬
ory, methodology, and verification of the fi¬
nite element numerical model AEXTIRC. The
model was developed to produce long numeri¬
cal simulations on the order of a year for very
large computational domains: for example,
the entire east coast of the United States. The
model was designed for high computational
efficiency and was tested extensively for both
hydrodynamic accuracy and numerical stabil¬
ity. Results of the tests are included in the
report.
AVAILABILITY OF REPORT: The report
is available through the Interlibrary Loan Ser¬
vice from the US Army Engineer Waterways
Experiment Station (WES) Library, telephone
number (601) 634-2355. National Technical
Information Service (NTIS) report numbers
may be requested from WES Librarians.
To purchase a copy of the report, call NTIS at
(703) 487-4780.
/.A PfocessKt Branch of the Coastal
at ^ US Anny EngiiM^r Waterways Eixperiment
Mr flie DRP wcHrk unit Dr, R, A. Luettich, Jr., Institute of
„.«rme Sciences. Uni^i;^M^?brt»Cteolina at 0^1 Hill, and Dr. J. I Westerink, Department
of CSvii Engineering, UiurertaQr bf Notre Dame, assisted Dr, Sclieffiio’ in developing the modeling
for fuid»»lnfoim E. ClaAMcNair,
I Ifi Manager* 0RF: at(^l) 6$4l207O, .
. — . ' . . ■■■■' • ' --
November 1992
Please reproduce this page locally, as needed.
REPORT DOCUMENTATION PAGE
form Approved
0MB No 0/04-0(88
PuWc fCportmg bufd'en iOf thu coU«ct(00 Ot mform^no^' n entmjied to ^ Hour o^( 'nporite. <nclodtAg !h« time lOr revif'^mg iniiruciiortt. ve^lfc^'nq eusting 04(« vOuk«.
gathering and maintaining the data needed, and <ompJetmg and reviewing the cohertion of information Ser>d <ommenw regarding th^ bijrden eMimaie or any other aspect of thi*
collectJOT) of information, including suggestions for reducing this burden to Washington Headquarters Services. Directorate for information Operations and Reports. UiS Jefferson
Davis Highway. Su»te 1204, Arlington. VA 222Q2-4302. and to the Office o* Managemer't and Budget. Paperwork Beduction Project {0704-0 1M). Washington. DC 20S03
1. AGENCY USE ONLY (Leave 6/anlc) 2. REPORT DATE
November 1992
4. TITLE AND SUBTITLE
ADCIRC: An Advanced Three-Dimensional Circulation Model
for Shelves, Coasts, and Estuaries; Report 1, Theory and
Methodology of ADCIRC-2DDI and ADCIRC-3DL
i. AUTHOR(S)
R. A. Luettich, Jr., J. J. Westerink,
Norman W. Scheffner
7. PERFORMING ORGANIZATION NAME(S) ANO AOORESS(ES)
See reverse.
3. REPORT TYPE ANO OATES COVERED
Report 1 of a series
5. FUNDING NUMBERS
WU No. 32A66
8. PERFORMING ORGANIZATION
REPORT NUMBER
Technical Report
DRP-92-6
9. SPONSORING /MONITORING AGENCY NAME(S) ANO AODRESS{ES)
US Army Corps of Engineers
Washington, DC 20314-1000
10. SPONSORING /MONITORING
AGENCY REPORT NUMBER
11. SUPPLEMENTARY NOTES
Available from National Technical Information Service, 5285 Port Royal Road,
Springfield, VA 22161.
12a. DISTRIBUTION /AVAILABILITY STATEMENT
Approved for public release; distribution is unlimited.
12b. DISTRIBUTION CODE
13. ABSTRACT (Maximum 200 words)
This report describes the theory, methodology, and verification of the
ADCIRC (ADvanced CIRCulation) finite element numerical models. ADCIRC-2DDI is
a two-dimensional, depth-integrated hydrodynamic model. ADCIRC-3DL is a three-
dimensional model that couples the external mode solution generated by ADCIRC-
2DDI to a locally one-dimensional internal mode solution. The numerical internal
mode solution can be obtained using either velocity or stress as the dependent
variable.
ADCIRC has been developed to simulate hydrodynamic circulation along
shelves, coasts, and within estuaries. To allow long numerical simulations
(on the order of a year) over very large computational domains (for example, the
entire east coast of the United States) , ADCIRC has been designed to have high
computational efficiency and has been tested extensively for both hydrodynamic
accuracy and numerical stability. The results of these tests are included in
this report.
14. SUBJECT TERMS
See reverse
15. NUMBER OF PAGES
141
Ifi. PRICE CODE
CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRACT
Of this PAGE OF ABSTRACT
UNCLASSIFIED UNCLASSIFIED
NSN 7540.01 -280-5500
Standard form 298 (Rev 2 89)
Pf»Hcnb<*q by ANSI Std 239 1%
298-102
7. (Continued)
University of North Carolina at Chapel Hill, Morehead City, NC 27514
University of Notre Dame, Notre Dame, IN 46556
USAE Waterways Experiment Station, Coastal Engineering Research Center,
3909 Halls Ferry Road, Vicksburg, MS 39180-6199
14. (Continued)
Bottom stress
Circulation model
Direct stress solution
Finite element method
Hydrodynamic model
Numerical model
Three-dimensional model
Two-dimensional model
PREFACE
The work described in this report was authorized and funded under Work Unit No. 32466,
"Numerical Simulation Techniques for Evaluating Long-Term Fate and Stability of Dredged
Material Disposed in Open Water," of Technical Area 1 (TAl), "Analysis of Dredged Material
Placed in Open Water," of the Dredging Research Program (DRP), sponsored by Headquarters,
US Army Coips of Engineers (HQUSACE). Messrs. Robert Campbell and John Lockhart, Jr.,
DRP Chief and TAl Technical Monitor from HQUSACE, respectively. Mr. E. Clark McNair,
Jr., Coastal Engineering Research Center (CERC), US Army Engineer Waterways Experiment
Station (WES), was DRP Program Manager (PM) and Dr. Lyndell Z. Hales, CERC, was Assistant
PM. Dr. Nicholas C. Kraus, Senior Scientist, Research Division (RD), CERC, was the Technical
Manager of the DRP TAl and Dr. Norman W. Scheffner, Coastal Processes Branch (CPB), RD,
CERC, is the Principal Investigator of Work Unit No. 32466. The numerical modeling goals,
concepts, and methodologies were developed by Drs. Norman W. Scheffner, Richard A. Luettich,
Jr., and Joannes J. Westerink. Development and in[q)lementation of the model were completed
by Drs. Luettich and Westerink.
This study was performed and the report prepared over the period of 19 July 1988 through
30 September 1990. Dr. Scheffner was under the administrative supervision of Dr. James R.
Houston, Director, CERC; Mr. Charles C. Calhoun, Jr., Assistant Director, CERC; Mr. H. Lee
Butler, Chief, RD, CERC; and Mr. Bruce A. Ebersole, Chief, CPB, RD, CERC. Ms. Janean C.
Shirley, Information Technology Laboratory, WES, edited the final report.
At the time of publication of this report. Director of WES was Dr. Robert W. Whalin.
Commander was COL Leonard G, Hassell, EN.
Additional information on this report can be obtained from Mr. E. Clark
McNair, Jr., DRP Program Manager, at (601)634-2070 or Dr. Notman W.
Scheffner, Principal Investigator, at (601)634-3220.
DTiC
iDl
Accesioo For
NTIS CRA&I
OTIC TAB
Unannounced
Justification
8y _ _
Distribution/
Availability Codes
Oist
I- .
Avail and/or
SoeciaC
1
CONTENTS
PREFACE . 1
LIST OF TABLES . 3
LIST OF FIGURES . 4
CONVERSION FACTORS, NON-SI TO SI (METRIC) UNITS
OF MEASUREMENT . 6
SUMMARY . 7
PART I: INTRODUCTION . 8
PART II: GOVERNING EQUATIONS . 11
Three-dimensional Equations for Nearly Horizontal Flow
in Cartesian Coordinates . 11
Three-dimensional Equations for Nearly Horizontal Flow
in a Coordinates . 13
Vertically Integrated, Two-dimensional Equations for
Nearly Horizontal Flow . 15
Mode Splitting . 20
Vertical Turbulent Closure . 22
PART III: EXTERNAL MODE SOLUTION . 24
Selection Considerations for the External Mode Solution . 24
Development of the Generalized Wave-Continuity Equation . 26
Bottom Stress Formulation . 29
Development of Weighted Residual Statements . 30
Time Discretization . 33
Spatial Discretization . 36
Solution Strategy . 48
Fourier Properties of the External Mode Solution . 50
Convergence Properties of the External Mode Solution . 56
Application of ADCIRC-2DDI to the English Channel
and Southern North Sea . 72
PART IV: INTERNAL MODE SOLUTION . 89
Definition and Applicability of a SDL Model . 89
Rationale for the DSS Technique . 90
Development and Testing of DSS Method No. 1 . 92
Development and Testing of DSS Method No. 2 . 115
Implementation of Wave^urrent Interaction in a DSS Model . 121
PART V: SUMMARY AND CONCLUSIONS . 130
REFERENCES . 132
2
LIST OF TABLES
No. Page
1 Constants for the Principal Tidal Constituents . 18
2 Summary of Error Measures for All Convergence Runs . 67
3 Am p 1 i tud e and Phase Relative to Greenwich for the Field Data,
Model Results . 81
4 Steady-State Bottom Stresses Computed Using Velocity
Expansions .
5 Analytical Solution for the Periodic Test Case . 105
3
LIST OF FIGURES
liOi Page
1 Modulus and phase of the propagation factor for PEFE and
PENSFD solutions . 52
2 Modulus and phase of the propagation factor for GWCEFE and
PESFD solutions . 54
3 Test problem geometry . 57
4 Grids used for the test problem . 58
5 Elevation and radial velocity components on the 6- by 8-
nodegrid with linear bathymet r y . 61
6 Convergence properties of the linear bathymetry test
cases . 68
7 Detail of the convergence properties of the linear
bathymetry test cases at low Cr . 69
8 Convergence properties of the quadratic bathymetry test
cases . 70
9 Detail of the convergence properties of the quadratic
bathymetry test cases at low Cr . 71
10 Fini te element grid and bathymetry used in the ADCIRC-
2DDI simulations of the English Channel and Southern
No rth Se a . 73
11 Time series of simulated and observed tidal elevations at
various coastal stations . 75
12 Time series of simulated and observed tidal speed and
direction at various stations . 78
13 Co-tidal charts for simulated constituents . 82
14 Drawings of the first eight Legendre polynomials in the
range -1 < <7 < 1 . 95
15 Vertical profiles of horizontal velocity for the steady
state test case for the VS and DSS^ . 101
16 Vertical profiles of the horizontal velocity magnitude for
the periodic test case for the VS and DSS^ . 107
17 Vertical profiles of the horizontal velocity phase for the
periodic test case for the VS and DSS‘ . 108
18 Ratio of the numerical and analytical bottom stress magni¬
tudes for the periodic test case for the VS and DSS ^ . 109
19 Ratio of the spectral to the analytical bottom stress phase
for the periodic test case for the VS and DSS^ . 110
20 Fourier series approximation to a square wave used to
represent an instantaneously applied wind for the
transient test case . 113
4
21 Time histories of numerical and analytical bottom stress
for the transient test case for the VS and DSS i . 114
22 Vertical profiles of horizontal velocity magnitude for
the periodic test case using the DSS| . 122
23 Vertical profiles of horizontal velocity phase for the
periodic test case using the DSS| . 123
24 Vertical profiles of horizontal velocity magnitude for the
periodic test case using the DSSg .
25 Vertical profiles of horizontal velocity phase for the
periodic test case using the DSSg . 125
26 Comparison of the DSS| and DSSg bottom stresses and the
analytical bottom stress for the pe r iodic test case . 126
27 Time history of bottom stress using the DSS| and the DSSg 127
for the transient test case .
5
CONVERSION FACTORS, NON-SI TO SI (METRIC)
UNITS OF MEASUREMENT
Non-SI units of measurement used in this report can be converted to SI (metric)
units as follows:
Multioh
degrees
feet
0.01745329
0.3048
To Obtain
radians
meters
6
SUMMARY
This report describes the theory, methodology, and verification of the finite
element numerical model ADCIRC, an ADvanced three-dimensional CIRCulation model
developed for the specific purpose of generating long time periods of hydrodynamic
circulation along shelves, coasts, and within estuaries. The intent of the model is to
produce long numerical simulations (on the order of a year) for very large computational
domains (for example the entire east coast of the US). Therefore, the model was
designed for high computational efficiency and was tested extensively for both
hydrodynamic accuracy and numerical stability. The results of these tests are included
in this report.
The ADCIRC model was developed by the Dredging Research Program (a) to
provide a means of generating a database of harmonic constituents for tidal elevation
and current at discrete locations along the east, west, and Gulf of Mexico coasts, and
(b) to utilize tropical and extratropical global boundary conditions to compute frequency
indexed storm surge hydrographs along the US coasts. The database of storm and tidal
surface elevation and current data is being developed to provide site-specific
hydrodynamic boundary conditions for use in analyzing the long-term stability of
existing or proposed dredge material disposal sites.
The overall intent of the DRP work unit is to provide a unified and systematic
methodology for investigating the dispersive or nondispersive characteristics of a disposal
site. These goals can be realized through the use of hydrodynamic, sediment transport,
and bathymetry change models. The ADCIRC model provides the tidal- and storm-
related hydrodynamic forcings necessary for site-specific site designation.
7
ADCIRC: AN ADVANCED THREE-DIMENSIONAL CIRCULATION MODEL
FOR SHELVES. COASTS. AND ESTUARIES
THEORY AND METHODOLOGY
PART I; INTRODUCTION
1. Interest in developing a more accurate technique for predicting sea surface
elevation and circulation in coastal areas has been spurred on by concerns relating to
navigation, shoreline flooding, pollutant transport, and sediment transport. A model
for computing the important features of circulation patterns driven by tides, wind,
atmospheric pressure gradients, and ocean currents must be broad in scope and size.
To simplify seaward boundary conditions, yet include important flow details, the model
must encompass large domains while providing a high degree of resolution in high-
gradient regions as well as in nearshore areas. This means that the model should
allow for the simultaneous solution of flow in continental shelf regions, coastal areas,
and in estuarine systems. The model should solve the three-dimensional conservation
equations [thereby resolving the vertical profile of horizontal velocity] instead of the
widely used depth-integrated conservation equations. This is necessary since it is
impossible to assume a relationship between bottom stress and depth-averaged velocity
that is generally valid for stratified flows, Ekman layers, and wind-driven circulation
in enclosed or semi-enclosed basins or in cases where wave orbital velocities or
suspended sediment concentration gradients are significant near the bottom.
f urthermore, it is impossible to assume values for momentum dispersion coefficients,
which are inherent in depth-integrated solutions, that are generally valid in complex
flows.
2. The requirements of very large domains, a high degree of horizontal
resolution in portions of the domain, and the resolution of i ’.pidly varying vertical
profiles of horizontal velocity place strenuous demands on even the largest
supercomputers. The goal in the development of ADCIRC (ADvanced three-
dimensional CIRCulation model) has been to bring together algorithms that are highly
flexible, accurate, and extremely efficient. These issues are closely interrelated and
have been emphasized in the selection of discretization techniques. The algorithms
that comprise ADCIRC allow for an effective minimization in the required number of
8
degrees of freedom for a desired level of accuracy, show good stability characteristics,
generate no spurious artificial modes, have minimal inherent artificial numerical
damping, efficiently separate the partial differential equations into small systems of
algebraic equations with time-independent matrices, and are capable of running months
to years of simulation while providing detailed intra-tidal computations.
3. The framework within which ADCIRC has been developed is a coupled
external mode - internal mode approach. This technique has proven to be successful
in past three-dimensional models and can significantly reduce the cost of three-
dimensional hydrostatic circulation computations. The governing equations and the
basic concept behind mode splitting are discussed in detail in Part II. The external
mode solution, which uses the well-known depth-integrated or shallow-water equations,
is discussed in Part III. Key features of the external mode solution include the use of
a generalized wave-continuity equation (GWCE) formulation and numerical
discretizations using the finite element (FE) method in space and the finite difference
(FD) method in time. Results are presented using the external mode solution as a
stand-alone, two-dimensional model on a quarter annular test case and the North
Sea/English Channel system. Part IV focuses on the internal mode solution'. During
the development of ADCIRC, a novel technique was discovered that replaces velocity
with shear stress as the dependent variable in the internal mode equations. The
resulting direct stress solution [DSS] allows physically realistic boundary layers to be
included explicitly in a three-dimensional model. This formulation of the internal
mode equations should be invaluable for modeling coastal and shelf circulation, in
which the bottom and surface boundary layers comprise a significant portion of the
water column, and for modeling processes that are critically dependent on boundary
layer physics such as wave-current interaction, sediment transport, oil spill moveuient,
ice floe movem energy dissipation, physical-biological couplings, etc. Thorough
descriptions of the DSS formulation and testing are presented in Part IV.
4. ADCIRC is being developed and implemented as a multi-level hierarchy of
models. A 2DDI (two-dimensional, depth-integrated) option solves only the depth-
integrated, external mode equations using parametric relationships for bottom friction
and momentum dispersion. A SDL (three-dimensional, local) option uses horizontally
decoupled internal mode equations to solve for the vertical profile of horizontal
velocity and to evaluate bottom friction and momentum dispersion terms for the
depth-integrated external mode solution. A 3DLB (three-dimensional, local,
baroclinic) option includes baroclinic terms as a diagnostic feature. Finally, the 3D
and 3DB options solve the complete internal mode equations for nonstratified and
stratified flows, respectively. At present ADCIRC-2DDI is fully implemented and
operational, ADCIRC-3DL is being tested, and other ADCIRC versions are under
development.
5. ADCIRC achieves a high level of simultaneous regional/local modeling,
accuracy, and efficiency. This performance is a consequence of the extreme grid
flexibility, the optimized governing equation formulations, and the numerical algorithms
used in ADCIRC. Together, these allow ADCIRC to run with order of magnitude
reductions in the number of degrees of freedom and the computational costs of many
presently existing circulation models.
10
PART II; GOVERNING EQUATIONS
Three-dimensional Equations for Nearly Horizontal Flow in Cartesian Coordinates
6. A survey of several recent review volumes (e.g., Heaps 1987; Nihoul and
Jamart 1987; ASCE 1988a, b; Davies 1989) indicates that the turbulent incompressible
Reynolds equations simplified using the Boussinesq approximation and the hydrostatic
pressure approximation generally form the basis for state-of-the-art numerical models
of coastal/shelf circulation. Although these equations describe fluid motion in three
dimensions, because of the simplification of the vertical momentum equations, they are
only correct for nearly horizontal flow (Koutitas 1987; Abbot 1990). Using a
right-handed Cartesian coordinate system these equations can be written as
du
dvi
W
dv
+ ^ + It “ ^
+ + V.
.dv
-f U^ + V
3y ^
dv
-dy
, dv
(1)
(2)
(3)
= - Pg
(4)
where
f = 2nsin^ = Coriolis parameter
g = acceleration of gravity
r = tide generating potential
1/ = molecular viscosity
p(x,y,z,t) = time-averaged pressure
/j(x,y,z,t) = density of water
Po = reference density of water
t = time
T = integration time scale for separating turbulent and time-averaged quantities
1
T
1/ ^
1
T
i;
r
u'u' dt - combined viscous and turbulent Reynolds stress
u'v' dt - combined viscous and turbulent Reynolds stress
11
T-zxC^.yi^.O =
r,y(x,y,z,t) =
7'yy(x,y,Z,t) =
'rzy(x.y»z.t) =
dw 1
^ ■ T
fT
I
0
T
u'w' dt - combined viscous and turbulent Reynolds stress
1
1 v'u' dt - combined viscous and turbulent Reynolds stress
dv 1
combined viscous and turbulent Reynolds stress
dw 1
u ~ T 1 “ combined viscous and turbulent Reynolds stress
(f> = degrees latitude
u(x,y,z,t), v(x,y,z,t), w(x,y,z,t) = time-averaged velocities in the x, y and z directions
u'(x,y,z,t), v'(x,y,z,t), w'(x,y,z,t) = departures of the instantaneous turbulent
velocities from the time-averaged velocities
X, y = horizontal coordinate directions
z = vertical coordinate direction
n = angular speed of the Earth (7.29212xl0‘® rad/s)
7. Using the vertical momentum equation, pressure can be eliminated as a
dependent variable from Equations 2 and 3, to give;
du , „5u , „5u . ^
, l/U C/U t
^ + v^ + - fv -
d
- r] + y^) - b. + m, (6)
^ + ''1^ + ^ I; + 8^ ■ ^] + ) - by + my (6)
where
bx
- S_ ^
by = S- I-
^ Po ^
(p-po) dz - baroclinic x - forcing
(p-Po) dz - baroclinic y - forcing
C(x,y>t) = free surface elevation relative to the geoid
Mtx = T-,
PoL
^ ^
35c by j
my = —
^ Po
- horizontal momentum diffusion
- horizontal momentum diffusion
Ps(x,y,t) - atmospheric pressure at the free surface
8. The solution of Equations 1, 5, and 6 requires the following boundary
conditions:
12
a. At the free surface,
Tzx — Tsx, Tzy — Tsy (J
where 7’sx(x,y,t) and rsy(x,y,t) are wind stresses applied at the water
surface.
b. At the bottom,
5h . „ 5h , 5hl
* = - iK + “ af + '
rz%/Po = Tix/Po = kUb, U,Ipo = 'Tby/po = kVb
(9)
(10a)
u = 0, V = 0
-h -I- Zo
(10b)
where Tbx(x,y,t) and Tby(x,y,t) are bottom stresses, Ub(x,y,t) and
Vb(x,y,t) are near bottom velocities, k is a slip coefficient and Zo is
the effective bottom roughness height (e.g., Zq = kg/30 where kg is the
physical bottom roughness). The physic^y correct no-slip condition,
Equation 10b, is often replaced by the slip condition. Equation 10a, to
avoid the need to numerically resolve the sharp vertical gradients of u
and V that exist near the bottom. A quadratic slip conmtion is
obtained by setting
k = Cd (ug+ vg)*^^ (11)
If the velocity profile is logarithmic between the elevation where Ub
and Vb are computed, (-h+zb), and the bottom, (-h+Zo), Cd can ^
defined rigorously as
Cd = {i ln[(zb-b)/(z,-h)l}-' (12)
where k is the von Karman constant. Often the quadratic slip
condition is replaced by a linear slip condition by setting k equal to a
constant.
c. At land boundaries normal flux is specified. Typically, this is zero for
a solid boundary or nonzero for a river boundary.
d. At open boundaries (either along the ocean or at rivers) the free
surface elevation, C(x>yit)» Is specified, a radiation boundary condition
is used to allow waves to enter and propagate out of the domain
(Davies and Fumes 1980; Rdd 1990), or the discharge is specified.
Three-dimensional Equations for Nearly Horizontal Flow in cr Coordinates
9. It is often useful to transform Equations 1, 5 and 6 into a bottom and
surface-following "a" coordinate system. By this means, numerical solutions of the
transformed equations maintain the same vertical resolution at each horizontal grid
point, regardless of variations in depth (Davies 1985; Blumberg and Mellor 1987). In
a general <r-coordinate system (where cr = a at the free surface and cr = b at the
13
bottom):
= X
Va = y
(T = a +
ta = t
where
H(x,y,t) = C + “ total water depth to the free surface
h(x,y) - bathymetric depth relative to the geoid
(13a)
(13b)
(13c)
(13d)
(The cr subscript is used to denote variables iu the new coordinate system.)
10. Derivatives are converted to the a-coordinate system using the chain rule;
d _ d , da a _ d (a-b) («<: . (<r-a) 8H 1 d
m - -d^w- m- ^
(14a)
d _ d , da d (a— b) , (<y-a) 5H 1 3
^ W H (a-b) ^
(14b)
d (a-b) d
(14c)
d _ d , da d _ d (a— b) f^C ■ (^^-a) 5H 1 5
M = m W = 'dt~ Tn ^
(14d)
11. The velocity component aligned in the a direction is defined as
do ^ (a^)
"" - 3t. TT
(a-b) 3r,
I ■ -B.* IS) g.)
12. The baroclinic forcings bx<,, by„, and the horizontal momentum diffusion
mxg, my^ in the a coordinate system become:
bxa - 2^^ %+ (p-Po) dor] + (cr-a) (16a)
L a
ha == (p-Po) + (o^-a) ^J^P-Po) (16b)
I <F
14
13. Substituting Equations 13 - 17 into Equations 1, and 4-6, and
rearranging terms gives the three-dimensional governing equations in the a-coordinate
system. Dropping the a subscripts for notational convenience, the transformed
equations are
5uH , dvE ,
(18)
dn . du ,
1
II
+ ^ - b. +
(19)
dv
+
(20)
(21)
14. The a equations use the same boundary conditions as the original
equations with the exception that w^= 0 at the free surface and at the bottom.
Vertically Integrated. Two-dimensional Equations for Nearly Horizontal Flow
15. The three-dimensional equations can be integrated over the vertical to
yield a set of two-dimensional equations for free surface displacement and
depth-averaged velocity. In conservative form these equations are;
aJH
+
5VH
W
(22)
auH , auuH . 5UVH
w~ + ~mr ~w~
- fVH = - H
+ g(C - a,)
15
(23)
+ M. + D. + B. + ISi - 32
Po Pq
m
5y
oy[po
+ My + Dy + By + ^ ^
^ ^ ^ Po Po
(24)
where
a = eflective Earth elasticity factor (a = 0.69)
Bx s - I bx dz - depth^ntegrated baroclinic forcing
-h
B
by dz - depth-integrated baroclinic forcing
-h
~ ~ ~ momentum dispersion
5Duv dPw
- momentum dispersion
Duu = [ uu dz, Duv = [ vu dz, Dw = [ vv
-h -h -h
r;(x,y,t) - Newtonian equilibrium tide potential
dz
Mx =
^ f
-h “h
dz - depth-integrated, horizontal momentum diffusion
n Q
My = ^ dz + ^ dz - depth-integrated, horizontal momentum diffusion
-h ^ -h
1
U(x,y,t) = w u dz - depth-averaged horizontal velocity
-h
1
V(x,y,t) = ^ V dz - depth-averaged horizontal velocity
u(x,y,z,t) = u - U - departure of horizontal velocity from depth-averaged velocity
v(x,y,z,t) = V - V - departure of horizontal velocity from depth-averaged velocity
16. In non-conservative form, the vertically integrated momentum conservation
equations are:
+ u|2 + vH - fV = - + g« - a,)]
+ h[«« + + b. + ^ - a.] (25)
16
(26)
It + + n’ = - + g(C - 0^)]
+ ^[m, + d, + b, + - ^]
17. The derivation of the Newtonian equilibrium tide potential, t), is presented
by Reid (1990). A practical expression for rj given by Reid is
1]{X,(p,t) = J Cjn fjn(to) hj(^) COS[2jr(t-to)/Tjn + jA + ^n(to)] (27)
n.j
where
Cjn = constant characterizing the amplitude of constituent n of species j (Table 1)
fjn(t) = time-dependent nodal factor
j = 0, 1, 2 - tidal species (j=0 declinational, j=l diurnal, j=2 semidiurnal)
Lo = 3 sin2(^) - 1
Li = sin(2^)
L2 = cos 2(0)
A, 0 = degrees of longitude and latitude, respectively
to = reference time
Tjn = constaiit characterizing the period of constituent n of species j (Table 1)
^n(t) - time-dependent astronomical argument
Values for fjn and can be computed from tables (e.g., Schureman 1941) or using
available harmonic analysis packages (e.g., Foreman 1977).
18. The gradient of ar) results in the effective tide-producing force. The
factor a accounts for the reduction in the field of gravity due to the existence of
small tidal deformations of the Earth’s surface called Earth tides. The value
a = 0.69 is the ratio of the theoretical period of the Earth’s wobble derived by Euler
(assuming the Earth to be a perfectly rigid sphere) to the observed period of the
Earth’s wobble (Reid 1990). (Therefore a is a global measure of the rigidity of the
Earth. For reference, a = 1 would correspond to a perfectly rigid sphere.) a = 0.69
has been used for modeling global ocean tides by investigators including Schwiderski
(1980) and Hendershott (1981).
19. Due to their computational efficiency, models based on the vertically
integrated equations have been widely used for modeling coastal, shelf, and even open
ocean circulation (e.g., Leendertse 1967; Wang and Connor 1975; Spaulding 1984;
Smith and Cheng 1987; Werner and Lynch 1987; Walters 1987; Vincent and Le
Provost 1988; Westerink, Stolzenbach, and Connor 1989; Signell 1989). All of the
physics contained in the original three-dimensional governing equations are embedded
17
Table 1
Constants for the Prindpal Tidal Constituents Cfrom Reid 19901
C
T
Species
Constituent
m
solar days or hrs* **
0
Mf
fortnightly lunar
0.041742
13.660791d
M„
monthly lunar
0.022026
27.554553d
Ssa
semiannual solar
0.019446
182.6211d
Sa
annual solar
♦*
365.2597d
1
Ki
luni-solar
0.141565
23.9344696h
0,
principal lunar
0.100514
25.8193417h
Pi
principal solar
0.046843
24.0658902h
Qi
elliptical lunar
0.019256
26.8683566h
2
Mj
principal lunar
0.242334
12.4206012h
S2
principal solar
0.112841
12.0000000h
N2
elliptical lunar
0.046398
12.6583482h
K2
luni-solar
0.030704
11.9672348h
*One lunar day = 1.035050 solar days or 24.8412 solar hours
**The annual solar tide is heavily dependent on seasonal heating and cooling of the
ocean, as well as radiation pressure.
18
in the vertically integrated equations if the bottom stress and the momentum
dispersion terms are specified correctly. Although more sophisticated approaches have
been developed for specialized conditions (Lynch and Officer 1985; Nihoul and Djenidi
1987; Tee 198’5’; Poon 1988; Jenter and Madsen 1989), bottom stress is usually
parameterized as a collinear function of the depth-averaged velocity, and momentum
dispersion is either neglected or represented as a "diffusion-like'' function of the
depth-averaged velocity (Bedford 1984).
20. Parameterized bottom stress relationships are typically quadratic in the
depth-averaged velocity and of the form
^ = Cf (U2 + U
Po
(28a)
^ = Cf (U2 + V2)i/2 Y
Po ' '
(28b)
where Cf is computed using one of the following relationships:
II
O
(29a)
C’
(29b)
C, =
(29c)
In Equation 29, is the Darcy-Weisbach friction factor, C is the Chezy friction
coefficient, and n is the Manning friction factor.
21. The depth-integrated lateral momentum diffusion terms are typically
lumped together with the momentum dispersion terms into a standard isotropic and
homogeneous eddy diffusion/dispersion model (Blumberg and Mellor 1987)
Mx + Dx =
My + Dy =
d^m
+
dy 2 dx. d y
pMD
hh j
ra^vH
dx 2
+ 2
d^VH
dy 2
a^UHl
(30a)
(30b)
where is a horizontal eddy diffusion/dispersion coefficient. Equation 30 is based
directly on a molecular diffusion analogy as applied to depth-integrated flow. Kolar
and Gray (1990) use a slightly simpler model that approximates Equation 30 as:
Mx + Dx
r^^uH
dy 2
(31a)
19
My + Dy =
t^MD
i^h2
r^^vH
~ur^
a^VHl
(31b)
where e“ j
is an eddy diffusion/dispersion coefficient that will generally not be equal
to E
MD
hr
22. For flows with horizontal length scales that are large compared to the
depth, Mx and My are negligible in the momentum balance in Equations 25 and 26
(Blumberg and Mellor 1987). D* and Dy are similarly small when the velocity profile
is nearly uniform over the vertical. In such flows Eh^ or Ehj are either set to zero
or kept at a relatively small value to provide stability to the numerical scheme. (The
latter must be done with considerable caution to ensure that the contributions of these
terms in the momentum equations remain small. Otherwise, the model solutions will
be artificially altered.) Conversely, when the velocity profile varies strongly over the
vertical, Dx and Dy may have a significant contribution to the momentum balance.
23. For tidad flows in relatively shallow, unstratified waters, depth-4ntegrated
computations that make use of the parameterizations given in Equations 28-31
appear to work reasonably well (although detailed studies of tidal constituent dynamics
indicate that all of the flow physics are not captured in two-dimensional simulations
due to the form of the bottom friction term (Westerink, Stolzenbach, and
Connor 1989)). However, in wind-driven flows, stratified flows, Ekman layers, or
when wave orbital velocities or suspended sediment gradients are significant near the
bottom, the simple parameterizations for bottom friction and momentum dispersion
given »above become entirely inadequate. Also, since the depth-averaged velocity may
be very different from the actual velocity at a specific elevation in th? water column
(particularly if flow reversal occurs over the depth), the use of the depth-averaged
velocity in a transport model (e.g., for sediment transport) may cause considerable
error in predicted transport patterns. Therefore, for many applications of practical
interest, a model based solely on the vertically integrated governing equations is not
adequate.
Mode Splitting
24. Unfortunately, numerical solutions of the three-dimensional governing
equations require substantially increased computer time and storage in comparison to
solutions of the vertically integrated equations. To help minimize this cost, most
three-dimensional models use some type of mode-splitting scheme. Mode splitting is
accomplished by solving the two-dimensional, vertically integrated, "external mode"
20
equations for the free surface displacement (and sometimes the depth-averaged
velocity). The external mode solutions are then used to force "internal mode"
equations that account for the vertical propagation of momentum. The internal mode
equations are solved for the vertical profile of velocity and the results used to
compute Tbx, Tby, Dx, and Dy for subsequent external mode calculations. Internal
mode equations have been generated by integrating the three-dimensional equations
over discrete layers in the vertical and then subtracting the equations for adjacent
layers (Simmons 1974; Sheng and Lick 1980), by subtracting the external mode
equations from the three-dimensional equations (Wang 1982; Sheng 1983; Davies 1985),
by differentiating the three-dimensional equations in the vertical direction (Tee 1979),
or by using the three-dimensional equations themselves (Blum berg and Mellor 1987;
Lynch and Werner 1991). (The internal mode equations and their solution are
discussed in detail in Part IV of this report.) Mode splitting allows the free surface
elevation to be evaluated with the computational efficiency of a vertically integrated
model. This can be quite important since the allowable time step for this
computation is often severely constrained by accuracy requirements or a Courant
stability criterion. Since the internal mode calculations are free from surface gravity
waves, the vertical profile of velocity can often be computed using a significantly
larger time step than the free surface elevation.
25. In effect, mode splitting replaces the parameterizations of bottom stress
and momentum dispersion used in a purely two-dimensional model with values
computed from the vertical profiles of velocity generated by the internal mode
equations. Therefore, the vertically integrated, external mode computations do not
require parameterizations of either bottom stress or momentum dispersion in terms of
the depth-averaged velocity. The only parameterizations maintained in the external
mode equations are for the horizontal momentum diffusion terms. These terms are
usually insignificant in the momentum balance, although for small-scale computations
horizontal momentum diffusion can be a physically important process. Most often the
horizontal momentum diffusion terms are retained only to provide numerical stability
and are parameterized with expressions identical to Equations 30 and 31, i.e..
52UH ^ a^VH'
ay 2 (hidy
(32a)
My
IjiM
Ehl
ra^vH
ax ^
+ 2
a^vH
~W^
+
a^uHi
WUy
or alternatively.
(32b)
Mx
ra^uH ^
a^uH'
ay 2
(33a)
21
(33b)
My = E^2
where E^j
ra^vH . a^VHl
and E^j are eddy coefficients for horizontal momentum diffusion.
Vertical Turbulent Closure
26. The internal mode equations require the parameterization of the vertical
turbulent momentum transport terms, and ray, (also called the vertical shear
stresses). These terms can dominate the momentum balance in portions of the
domain and it is therefore critical to find an adequate closure scheme. Turbulent
closure has been and continues to be the subject of considerable research. Recent
summaries of this work include Mellor and Yamada (1982); Rodi (1984, 1987);
Ferziger (1987); Johns and Oguz (1987); and ASCE ( 1988a, b). The most general
approach is to solve transport equations for the turbulent velocity correlations that
make up the turbulent stresses (stress/flux models). However, this adds considerably
to the computational burden of a three-dimensional model. Models based on this
technique have had little testing and virtually no application to geophysical flows
(ASCE 1988b). Also, it appears that these models offer no decisive advantage in
shear flows (Launder 1984). Alternatively, the vertical shear stresses can be
parameterized in terms of the mean velocity field using eddy viscosity relationships of
the form
Tzx _ p ^
(34a)
= Ey
po ^ m
(34b)
On dimensional grounds the vertical eddy viscosity Ey should be proportional to a
velocity scale v multiplied by a length scale I, both of which are characteristic of the
turbulent motion. Particularly simple expressions such as the Prandtl mixing length
model can be found for v and I for boundary-layer type flows (Rodi 1987). In more
complex flows, v has been related to the square root of the total turbulent kinetic
energy, k. The terms k and / (or some combination of k and I such as can
be solved for using quasi-empirical transport equations or specified using empirical
algebraic expressions. The primary limitations to the eddy viscosity approach are its
inability to simulate counter gradient transport or to a^'count for nonisotropic
turbulence. A third choice for expressing the turbulent stresses lies between the
stress/flux models and the eddy viscosity models in complexity and potential for
22
representing co>^.Dlex flows. In this approach algebraic expressions (approximations to
the transport equations used in stress/flux models) relate the vertical stresses to k and
I (or e) without the use of an eddy viscosity hypothesis.
27. Eddy viscosity models are by far the most widely used method for
representing vertical momentum transport in coastal flows. These models can be
expected to work reasonably well in such applications, since the water column is
typically dominated by the bottom and surface boundary layers.
23
PART III: EXTERNAL MODE SOLUTION
Selection Considerations for the External Mode Solntion
28. A basic objective in the development of ADCIRC is to provide the abibty
to perform computations on very large domains. This requires selecting algorithms
that satisfy interrelated requirements of a high level of grid flexibility, accuracy, and
efficiency. To ensure a high degree of solution accuracy, the discretization scheme
must have numerical amplitude and phase propagation characteristics that are nearly
identical to the analytical characteristics even for relatively poorly resolved
wavelengths (e.g., good correspondence down to at least A/ Ax = 20, where A is the
wavelength and Ax the grid spacing). Furthermore, solution accuracy requires that all
wavelengths with significant energy, (e.g., as generated in regions of rapidly varying
flow, geometry, and/or topography), be well-resolved. A high degree of solution
efficiency requires that the algorithm minimizes both the number of degrees of freedom
and the operations required per degree of fireedom per time step. Minimization of the
number of degrees of freedom is constrained by the need to provide resolution on a
localized basis and is highly dependent on the accuracy and the grid flexibibty of the
numerical scheme.
29. Because grid flexibility is pivotal to solution accuracy and efficiency,
various strategies have been devised to allow variations in grid size over a model
domain, A nested grid approach offers one solution. However, unless the grids are
coupled, this approach cannot properly account for flow interactions between the
various grids. Stretched FD grids offer the possibility of providing local refinement
within a single grid. However, cell aspect ratio requirements limit the degree of grid
size variability. Furthermore, since cell size in the x direction is fixed for all y
locations for a given x and vice versa, portions of the domain are often over-refined.
Boundary-fitted FD schemes that utilize conformal mapping techniques allow the land
boundaries to be well-represented in addition to offering local refinement possibilities.
However, these techniques suffer from the same shortcomings as stretched FD
approximations and often significant difficulties are encountered in finding a suitable
transformation function for complex geographic regions. The FE algorithms based on
triangular elements are highly flexible and can provide local refinement in a systematic
and optimal fashion. In fact, circulation computations for tides and storm surge in
the Gulf of Mexico (Westerink et al., in press) have been achieved with cell area
ratios greater than 1 to 15,000.
24
30. Algorithm accuracy per degree of freedom is another critical issue in the
selection of an external mode solution algorithm. FD schemes were successful fairly
early in their development owing to the use of the staggered or C grid approach
(Hansen 1956; Leendertse 1967). Early FE schemes were plagued with severe spurious
modes that required the heavy-handed addition of non-physical dissipation and
resulted in very poor accuracy characteristics (Gray 1982). It was not until the
introduction of the wave-continuity equation (WCE) formulation that robust and
highly accurate FE schemes emerged (Lynch and Gray 1979). The WCE formulation
is based on the rearrangement of the continuum equations prior to any spatial
discretization. Extensive numerical testing has demonstrated that FE-based WCE
solutions produce very accurate results (Lynch and Gray 1979; Lynch 1981; Walters
and Carey 1983; Walters 1983 and 1984). It has also been shown that the
fundamental success of the WCE FE scheme lies in its ability to propagate 2 Ax
waves (Platzman 1981; Foreman 1983). (This is also the reason why the C grid FD
solutions are successful.)
31. Finally, overall algorithm efficiency is essential in the selection of an
external mode solution. In general, implicit methods are more useful in long wave
computations than explicit methods, particularly when small cells or elements are used.
However, the use of implicit methods typically results in time-dependent matrices that
must be reassembled and re-solved at every time step. This increases the
computational burden significantly. The FD methods overcome this problem by
implementing an alternating direction implicit (ADI) type approach that reduces a
two-dimensional problem to a sequence of one-dimensional problems, resulting in
significant computational savings for large problems. It is not possible to apply the
ADI approach to FE-based methods. However, a WCE FE-based solution has been
formulated that decouples the solutions for elevation and velocity and allows the use
of time-independent matrices for the elevation solution and diagonal matrices for the
velocity solution. These features have produced a highly efficient WCE FE solution
called the generalized wave-continuity equation (GWCE) formulation (Kinmark 1985).
32. Careful consideration of the requirements for grid flexibility and a high
level of accuracy and efficiency led to the selection of the FE-based GWCE
formulation for the external mode solution in ADCIRC. Extensive analysis, testing,
and field applications of the GWCE during the past decade have demonstrated the
unparalleled capabilities and robustness of the scheme.
25
Development of the Generalized Wave-Continuitv Equation
33. The GWCE formulation is a specifically designed WCE formulation that
yields a discrete system of equations with time-nndependent matrices. Time-
independent system matrices are critical in minimizing the computational cost for
finite-element-based solutions due to the expense of both the matrix assembly and
decomposition steps. The GWCE is based on the primitive depth-integrated
continuity equation, Equation 22, and the primitive depth-integrated conservation of
momentum equations in conservative form. Equations 23 and 24. The primitive
continuity equation is differentiated vrith respect to time to yield:
^ ^t^ wdj - ^
The primitive momentum equations are differentiated with respect to x and y,
respectively, and rearranged as:
aJUH WVH
•f Mx + Dx + Bx + ^
Po Po^
auvH avvH
+ M, + D, + B, +
Equations 36 and 37 are then substituted into Equation 35:
d I gUUH gUVH
^ * dx 9y
+ fVH - H (^ + g« - av)\ + M,
+ D, + B, + ^
Po
- H ^ + g(C - 07?)] + My + Dy + By + ^ - = 0 (38
Finally, the primitive continuity equation is multiplied by a constant, Tq, and added
to Equation 38:
Thx . Tttxl . d , WVH ^VH rrrtT
+ 7-oUH} + 7^ {- - - fUH
26
- H ^ 1^ + g« - <«,)) + M, + D, + By + + ToVH} = 0 (39)
34. The advective terms in Equation 39 are in conservative form. Our
experience indicates that if these terms are put into non-conservative form, improved
numerical stability is obtained when advection is dominant in the global or local force
balance. The advective terms in the GWCE are reformulated by expanding the
derivatives and substituting in the primitive continuity equation, Equation 22.
+ -. ^ + If {u - UH II - VH |S + fVH - H ^ + g(C - a,)l
+ M. + D. + B. + ^ - + r.UH} + I (V § - UH H
- VH || - fUH - H ^ [fe + g(C - c,„)] + My + Dy + By
+ ^ ^
35. The lateral closure model in ADCIRC is the simplified eddy viscosity
model of Kolar and Gray (1990), Equation 33. Substituting this into Equation 40
gives:
+ -.|^ + If {U ^ - UH H - VH H + fVH - H I; [^ + g« - a,)l
+ Dy + B, + + t„UH} + {V § - UH H - VH || - fUH
- H ^ + g« - 09)1 + Dy + By + 2^ - + ^.VH}
. d ,p /d^UH , a’UH,, . 3 ,p ,3^VB
where Ej^j is the generalized lateral diffusion/dispersion coefficient. For the 2DDI
option, E}j2 represents the combined effects of both lateral diffusion and dispersion.
Therefore Ejjj = and D* and Dy are both set to zero. For the three-dimensional
ADCIRC options, Ej^j represents only lateral diffusion. In these cases, Ejjj = E^j,
and Dx and Dy are explicitly computed from the internal mode solution. It is
assumed that E^^j is constant in time and space and that it has a value of zero on
the boundaries of the domain.
36. The lateral diffusive/dispersive terms in Equation 41 can be conveniently
rearranged to decrease the functional continuity requirements for the symmetrical weak
27
weighted residual formulation from back to C® as is the case for the GWCE
formulation without any lateral closure model (Kolar and Gray 1990). Rearranging
the spatial derivatives of the lateral diffusive/dispersive terms in Equation 41 gives:
0 It ■ H ~ I7 + ® H
+ D, + B. + -h roUH} + I {V 1^ - UH 1^ - VH 1^ - fUH
- H ^ [Es + g(C - at,)] + Dy + By + I^-I^ + toVH)
. ,, td^ fdVE , 5VHm . ^ td^ ,5UH . 5VHm _ ^
+ Ehj (^3^ + ^)] + Ehj + ^)] = 0 (42)
The primitive continuity equation, Equation 22, can be used to substitute for the
divergence of flux in the lateral diffusion/dispersion terms in Equation 42 to give:
+ ^.^ + If {u § - UH |S - VH ^ + fVH - H 1^ + 6(C - a,))
- & + Dx + B, + + r,UH}
+ ^ {V 1^ - UH - VH - fUH - H ^ + g« - a,)l
37. Equation 43 can be solved in conjunction with the primitive conservation
of momentum equations in either conservative or non-conservative form. ADCIRC
uses the non-conservative momentum equations, Equations 25 and 26. Incorporating
the same simplified eddy viscosity model into the non-conservative momentum
equations gives:
28
38. The bottom stress in Equations 43 - 45 is expressed using a drag tensor
similar to that proposed by Jenter and Madsen (1989):
and 7 is the angle measured counter clockwise &om the depth-averaged velocity vector
to the bottom stress vector.
39. Defining
P= f + r*sin(7) (47a)
ri = r*cos(7) (47b)
aiiH substituting Equations 46 and 47 into Equations 43 — 45 gives the GWCE and
momentum equations in final form:
+ -«^ + H ^ ™ li - VH f + f'VH - H (^ + g« - a,))
- Eh. & + D, + B. + + (r<,-ri)UH}
+ I? H f - f'UH - H I + g« - a,))
- Eh. + Ey + B, + 2^ + (r,-ri)VH} = 0 (48)
- t'V = - If [^ + g(C - on)]
H Eh.
ra^uH . a’uH
ru-cin , u vni , v/x . . 'sx
l^TT- + ^T-1 + H- + BT + KH
It + U I? + V + f'U = - If (^ + g« - o.,)!
+ H Eh. l|;7® + |r^l + ^
40. In the 2DDI option, the bottom stress and depth-averaged velocity are
assumed to be co-4inear (7 = 0). Cf is specified directly as an input parameter or
computed using one of the relationships given in Equation 29. In the three-
dimensional ADCIRC options, 7 and Cf are computed using rb* aJ^d rby from the
internal mode solution. As noted above, 7 is the angle measured counterclockwise
from the depth-averaged velocity vector to the bottom stress vector. Cf is determined
as:
Po(U‘ + V’)
(51)
It is easily shown that Equations 46, 47, and 51 introduce the bottom stresses
computed in the internal mode solution directly into the external mode equations.
Development of Weighted Residual Statement
41. To develop a Galerkin weighted residual statement for the GWCE, Equation
48 is weighted by the interpolating basis function, and spatially integrated over
the interior domain, D, giving:
^ -N (52)
where
<o,6>n = ^/ 0 6 dn
n = the global domain
N = number of nodes in the spatial discretization
A. = U If - UH |2 - VH ^ + f'VH - H (^ + g« - a,)]
- Ek, 1^ + D, + B, + + (r<,-Ti)UH (53a)
A, = V If - UH 1^ - VH f - f'UH - H (Ei + g(C - a,)]
- Ehl + ^ + (ro-rt)VH (53b)
Applying Gauss’s theorem to the integrals in Equation 52 that contain spatial
derivatives gives:
<|t^> ~ - <Ay, |^>n
= - j, [Ax«nx + Mnyl^i " N (54)
30
where F is the boundary of the domain fi. The direction cosines are defined as:
Onx = cos{^x) (55a)
Ofay = cos(^y) (55b)
where and &y are the (spatially varying) angles measured to the outward normal at
any point along the boundary from the positive x and y axes, respectively.
42. Using the conservative form of the momentum equations, Equations 23 and
24, recasting the advective terms in Equation 53 into conservative form, and using the
simplified lateral diffusion model, Equation 33, Ax and Ay can be written as:
A, =
5UH
- ^h2 +
Ay - - l^ha +
d^VH.
'dp-^ ~
+ ^oUH
+ r,VE
(56a)
(56b)
Substituting Equation 56 into the line integral in Equation 54 and assuming that Ejjj
is zero on the boundary. Equation 54 becomes:
It “ lx ~ ly^
“ ® S 8(C - <w?)l - Et, ^ + D, + Bx + ^>n
- <V ^ - UH - VH f - f'UH - H I; (fe + 6« - a,)] - E„ ^
+ Dy + By + + (r<,-ri)VH. ^>a = - I ||{UHa„x + VHa„,)
+ To(UHOfnx + VHttny)]^! dF i = 1, ...N (57)
43. The terms that involve partial derivatives of the barometric pressure,
surface elevation, and Newtonian equilibrium tidal potential can be written as:
“ If + 8(< - Of)] = Sk ^ + I ^ + 6H (^ - OT,) (58a)
H I - “^)1 = Si I + § I + SH I (^ - «5) (58b)
44. The normal flux across the boundary is defined as:
Qn = UHOnx "b VHOfny (59)
45. The line integral in Equation 57 is non-zero only on flux-specified
boundaries, F^. Using the specified normal flux Qn* for Qn, and substituting
Equations 58 and 59 into Equation 57 gives the final symmetrical weak weighted
31
residual statement for the GWCE:
+ <To^,<t>i>n + <g^^> ^>n + ^>n + ^>n
+ Eha<^, ^>n = <u|^. ^>n + <v|^. ^>n + <Wx, ^>n
+ <Wy, ^>n - ^ (#"* + ^oQn*) dr i = 1, ...N I
where
w, ^ - UH H - VH 1^ + f' VH - I g! - gH (-^ - a,) + D, + B.
+ ^ + (ro-ri)UH
(6U)
W, . - UH W - VH 1^ - f'UH - § gH - an) + D, + B,
+ ^ + (ro-rJ)VH
(61b)
46. The weigh.ed residual form of the conservation of momentum equations is
obtained by weighung Equations 49 and 50 by and integrating over the domain Cl:
<It + « W + V f + t'U + I [?^ + g(C - ar,)]
- H Em - ^ + riV, ^^>a = 0 (63)
Applying Gauss’s theorem to the lateral diffusive/dispersive terms in Equations 62 and
63, and recalling that Ej^, equals zero on the boundary, gives the symmetrical weak
weighted residual form of the momentum equations:
32
(65)
- I?; + 6(C - <«/)l - ^ + riV, #,>„
“ + V ^i>ii + ^i>n
47. The GWCE equation is discretized in time using a variably weighted,
three-time-4evel, implicit scheme for the linear terms (i.e., those terms on the left side
of Equation 60). Wx and Wy are treated explicitly. The time derivatives that
appear on the right side of Equation 60 are evaluated at two known time levels. The
time-discretized GWCE is:
_ 2 A 4. /k-l /k*l _ A-l
— + To<^ . , ^^i>n
+ «i [<g^^ . ^>n + . ^>nl
+ «2 [<g^^ > ^>n + <g^^ > ^>nl
+ “3 l<g^^ . ^>n + <g^^ • ^>nl
= . <V^(^), f >„
+ <w|, ^>n + <wf,
+ ToOr^Qn^E + T'otksQn**) dP i = 1, ...N
where
At = time step
k+1, k, k-1 = future, present, and past time levels
ttj, Oj, Oj = time weighting factors
The time weighting factors are selected so that:
OTj + aj + Oj = 1
o, = 03
Rearranging Equation 66 gives:
33
(1 + ^1>„
+ a,gAt=[<h^ , ^>„ + <h^ , ^>(J
+ -a^M [<^-, ^ <|-,
= 2 <<‘, h>a + - 1) <C'“'. ^i>n
- ajgAl2[<h^ , ^>a + <h^ , ^>n)
- . ^>n + <1>|^ . ^>nl
+ l<r“. ^>« - <t'-
+ At<Ui‘«i' - ^>„ + AKVHi" - ^>a
+ At’<W!S, ^>„ + At=<w!;, ^>n - At’F, i
i=l, ...N
where
= I + ToaiQli;! + ToazQS* + ToChQU^) (t>i dr (68]
48. The symmetrical weak weighted residual form of the momentum equations
are discretized in time using a two-time-level implicit Crank-Nicolson approximation
for all terms except the diffusive terms, which are treated with a variably weighted,
two-time-level implicit scheme and the advective, dispersive, and baioclinic terms,
which are treated explicitly:
+ Z <n'‘(U'‘*' + U''), 0i>n - + V-), *,>„
+ E„ [ft <42hi-',
, ft , ft
= - z <1:1?!^ + - (jS)-'’ '^On
- Z 4 [?! + sfC- - - (^h)'. *>n
^TTk 3U'‘ . vk 5U‘‘ 2 ^ Bil , ^
_ <U ^ + V ^ , ^i>n + + IJJ, <f)i>^
i=l, ...N
< — At + Z + yk), ft>„ +<-|l'‘{u'“' + u^), ^,>n
34
+ Eh! [^1
, A <^"; |(^>„ . /*!
= - J <|^?!^‘+ ece- - oV“‘)l - (gi)-', h>«
- 2 <1^ 1?! + E«‘ - “-^ll - ^»n
- <u'‘ + Vk *,>„ + <g{ + B|, ^,>„ i=i, ...N (70)
where /3i and 02 are time-weighting factors at the future and present time levels.
These factors are selected so that
01 + 07 = ^
Rearranging Equations 69 and 70 gives:
<(1 + firjk) Uk-J <ln>a - <f'kvk‘;
= <(1 - ^rjk) Uk, #,>„ + <£'kvk, ^i>„
- ¥ <® 6«'‘" - - (^)'‘"> '^»n
- At <Uk |2 + ly ’ 't'i>a ‘=k> ■N
<(1 + ^rjk) yk'l ^i>„ + ^ <f'kuk>i, 0i>„
= <(1 - lirjk) yk, 4,i>„ - ^ <t'kuk,
- ^!Eh!A‘ !<^‘. |(fe>0l
- ¥ <|f !?f‘+ eK"*' - “V"‘)l - (gt)'"'. h>a
-¥<ht^
35
i=l, ...N (72)
- At <U'‘ + ]^» ^i>n
49. There are two differences in the solution of Equations 67, 71, and 72 in
the 2DDI option and the three-dimensional options. In the 2DDI option, the friction
parameter (Cf or one of the parameters in Equation 29) is specified in the model
input. The dispersive terms are included with the lateral diffusive terms by
eliminating Dx and Dy from Equations 67, 71, and 72 and setting E2h = Eah- In the
three-dimensional options, Cf, 7, Dx, and Dy are computed from the most recent
internal mode solution. In flows where the velocity reverses direction over the depth,
it is possible for the depth-averaged velocity to be zero while the bottom stress is
nonzero. In this case the drag coefficient computed in Equation 51 becomes infinite.
To prevent the numerical difficulties that this causes, an upper limit is set on the
computed drag coefficient. If this limit is exceeded, 7 and Cf are set to zero and the
bottom stress computed in the most recent internal mode solution is passed directly to
the external mode equations. In the GWCE, — and ^ determined in the internal
Po Po
mode solution are subtracted from W* and Wy, respectively. In the momentum
equations, and are subtracted from the right-hand side of the
corresponding equation. This modifies the final terms in Equations 71 and 72 to
+ HE respectively.
Spatial Discretization
50. In order to complete the conversion of the governing partial differential
equations into systems of algebraic equations, the FE method is applied to the time-
discretized form of the symmetrical weak weighted residuwl equations developed in the
previous section. Specifically, elemental approximations to the variables are
substituted into Equations 67, 71, and 72, the elemental equations are summed over
the global domain, and the required degree of inter-element functional continuity is
enforced. Interpolating basis with at least C® functional continuity are required to
discretize most of the dependent variables. Departures from this are noted below.
51. In all linear terms, surface elevation, velocities, and depth are
approximated over each element as:
36
(73a)
<‘‘5 E <t'i
j * 1 ■' •*
^ E Uf
j * i
V'‘ 5 E
j =1
J1 e 1
(73b)
(73c)
h 2! E h: (73d)
j . 1 J *
where equals the number of nodes per element. In nonlinear terms and certain
forcing terms, the entire term may be interpolated over the element as described
below.
52. The nonlinear and forcing terms in the GWCE are approximated as
follows.
a. The Coriolis parameter and the fluxes in the Coriolis term are
appro.'dmated by:
f'’‘(UH)>‘ ~ ”e^ (f'UH)f 4>-, (74a)
i =1 * *
f/k(vH)k !v (f'VH)^ 4>.
j « 1
fe. The finite amplitude component of the free surface gradient is
approximated by:
«7)k ? “e (OJ
j * 1
(74b)
(75)
c The combined barometric pressure and Newtonian tidal potential term
is approximated by:
(5?i -
Vog " - j M Vog
The total depth factor in this term is evaluated using an L:
approximation:
2; ^ “e^ H5
- nei j J
d- The surface stress terms are approximated b3';
Ue 1
(Isi)k ^ 6.
Vo - j =1 Vo
^77)
(78a)
(78b)
37
e. The bottom stress and Tq terms are approximated by;
(rik- ! “e’ [(rj - ro)HUlf A
j * 1 ^
(79a)
(Ti'‘- To)H'‘V'‘ 5 E l(r; - r„)HV)} (79b)
j * I
In the three-dimensional model options, if the computed friction
coefficient exceeds the maximum allowable value, the bottom stress
terms are approximated directly by:
k Del
Iki ~ JJ {'Ibs)k A.
Po = jTi Vo
N *^1;^ (Ibi)k A
Po = jTi Vo h
(80a)
(80b)
f. The dispersive terms are broken up into their components Duu, Duv
and Dw, (defined in Part II), and discretized as;
S “e_ ^
I^Dvv)'- S ^
(81a)
(81b)
(81c)
(81d)
g. The baroclinic terms are not included in either ADCIRC-2DDI or
ADCIRC-3DL. Therefore the discretization of these terms is not
considered here.
h. The velocities that multiply the time derivative components of the
non-conservative advective terms are approximated using L2
interpolation:
^ = -
1 “el ,
el j = t •'
(82a)
v‘ = E 5^ ”e VJ (82b)
The free surface elevations that appear in these terms are
approximated using the standard Co approximation. Equation 73a.
The spatially differentiated components of the non-conservative
advective terms are approximated by
mi , , ®el .
(UH §)‘ = (UH)J, ,E U} ^
(83a)
38
(VH 1^)“ ! (VH)5, U} ^ (83b)
(UH 5 (UH)5i ”s‘ V8 ^ (83c)
(VH |^)‘ S (VH)J, vy ^ (838)
where
f ^ 6 1
(UH)Si H ^ E (UH)} (84a)
“el j = I
1 Q e 1
(VH)5, e i .E (VH)8 (84b)
53. Substituting the approximations in Equations 73 - 84 into Equation 67
and summing over the elements gives:
He 1
H
E
e 1 »l
+
+
j?. [d +
<.igAl2[< ||i, ^>0,,+ <^S ^>nj
^ ^>0., + «i’‘
- “i8^d[<^E^ h„0„<y ^>n„+ j M.<i
- h.0.<S- <”s_
^ «4- ^>0, + «S- f. ^>«J
At l<US.(fl - <f-) 4>j, + <V‘,«5 - <?-') fi>„J
- Att[<(UH)5.u5 gi>„^,+ <(VH)S.US
+ <(UH)S.VS[ + <(VH)‘,V}
+ Attl<(f'VH)8 - <(rUH)f ^>„J
- ^ 1<(05 ^>8., + <«’)? f >»,.i
- «AttHj.l<(^ - a,)) ^ >„,,+ <(^ - m,)) |i. fi>„.,l
+
+
39
+
- At'{<l(rj - r.)HU|5 + <[(rj - r.)HV]; 0^,
where
M = the total number of elements
flgj = the elemental domain for the element, el.
Equation 85 can be rewritten as;
M
E
e 1 *t
nei
E
j *1
[(1 +
= "e‘ [2m(|)<J + (2^ - - a,gAt^M(?)<J
+ IM^SlCCf - Cf") + - <}-)j
- At>(M(t>(UH)‘.uf + M(f(VH)5,u} + m!?(uh)5,v} + m(?)(vh)5,v;i
+ At*(M(p(f'VH)J - M(p(f'UH)Jl - MjP(Oj
- At*{M(T)[(r; - r.jHU]} + M(f)[(rJ - rojHV]})
- At>[M(pD„„J + M(f)D„vf + M(pD„,J + M(pDv,{]]
- At’F,
i=l, ...N
where
Del
+ < E
m * 1
Mi?) ^ KK t-
Mi?) s
(86)
(87a)
(87b)
(87c)
40
mW . <^. ^>0^,
mIP ^ <^i> ^>0.,
mSP = <^j.
Note that the elemental matrices, Mi}) , mW, m}}) . Mif,) and Mif) are
symmetrical and that M^p, and M^p are non-symmetrical.
54. The fully discretized GWCE can be written in a compact form as:
E {Cf 1} = E i=l, ...N
c 1 “1 j « 1 ■* e 1 »1
(87d)
(87e)
(87g)
(87h)
where
Mff®® s (1 + ^)Mi}) + a^At>Mif) + 56^ m}}) (
pewcE ^ ”1;'^ [sMi})(} + - l)Mi})<}-' - a,gAt’Mif)<}
- arfAt=Mif)<}-‘ + Mif)<}-‘
+ At (Mi})uk,(C} - <}■•) + Mif)vk,«} - <}-))
- At^Mif)(UH)‘,U; + Mif)(VH)iiU} + Mif)(UH)5,V5 + Mif)(VH)5iV5l
+ At*IMi})(f'VH); - Mif)(f'UH);] - Mif)(OJ
- 8A‘’Mif)(^ - + At=[MiP(^)} + Mif)(I“)}l
- AtHMi})[(n - t„)HU)} + mWKtJ - r,)HVl}}
- At’(Mi})D„„5 + Mif)D„v} + Mif)D„,5 + Mif)D„}l] - At>Fi
i=l, ...N (
In the three-dimensional options, if the computed friction coefficient exceeds the
CWCE
maudmum allowable value, the friction term in the right side load vector Pj is
slightly different:
41
pGWCE ^ “g‘ + (2^ - - o^At“M(?)c5
+ At [MipUk,((:V - <lf->) + M(®)v5,«j[ - C|f')l
- At’[M(t)(UH)‘,Uj[ + Mlf)(VH)‘,U} + M®(UH)‘,V5 + M(^)(VH)k,V}l
+ At’IM^JWjf - Ml?){roH)J) - mSJ) (05
- gAtm)(^ - <»,)5Hk, + At=lM(J)(lM)k + m(?)(2«)})
- At“{M(J)[(^)5 - r,(UH)5] + m(®)[(J^)5 - t„(VH)51}
- At>(Mjt'D„„} + m(®)d„,5 + mWd„,5 + m(^)Dv,5]] - At»Fi
i=l, ...N (91)
55. Global assembly and enforcement of the functional continuity
requirement leads to the following global system of equations:
E [ ] { 8Ci*‘ } = { } j=l. •••N (92)
j M ■* J
where
g^GWCE _ global banded system matrix
gpGWCE _ global load vector
= the global nodal elevation vector
56. The fully discrete form of the momentum equations is obtained from the
time-discretized symmetrical, weak weighted residual form of the momentum equations,
Equations 71 and 72, as follows.
g,. The local acceleration terms are interpolated using Equations 73b and
73c.
The friction terms are approximated by:
j = i
(93a)
j * 1 •'
(93b)
42
(93c)
Del
j • i ■* ■'
(93d)
In the three-dimensional model options, if the computed firiction
a)efGcient exceeds the maximum allowable value, the bottom stress
terms introduced on the right-hand side of the momentum equations
are approximated by:
ttel
(94a)
(94b)
c. The Coriolis terms are approximated by:
^ € 1
f.kukM ^ 2 f/kuk*i ^
j * 1 '* ■' •*
nel
f/kvk»i ~
j = 1 ^ ^ ^
nei
f/kyk 2; £ f/kuk ^
j = 1 •'
nei
f/kyk N S f/kyk ^
i «i ■'
(95a)
(95b)
(95c)
(95d)
d. The lateral diffusive/dispersive terms are approximated using
Equations 73b and 73c for velocity and Equation 77 for total depth.
e. The barometric pressure and Newtonian tidal potential are
interpolated using Equation 76.
f. The surface elevation is approximated using Equation 73a.
g. The surface stresses are evaluated as:
n e 1
/'^SX'ik*! ni V
W
(96a)
n e 1
^ J?,
(96b)
n e 1
(96c)
^ € 2
(96d)
h. The advective terms are approximated by:
43
r = V,V£; U5 ^ (97b)
u^ t USi^E Vf ^ (97C)
V"' W ^ VS, “2, Vj ^ (97d)
where U^i and V^j are defined in Equation 82.
i. The dispersive terms are broken up into their components Duu, Duvi
and Dw, (defined in Part II), and discretized as:
(h f “"“f ^
(h f ” (hIi)''”?. “-i ^ (»*'>)
(h f®*')
(h f s (0^:)““?,
where Hei is defined in Equation 77.
j. The baroclinic terms are not included in either ADCIRC-2DDI or
ADCIRC-3DL. Therefore the discretization of these terms is not
considered here.
57. Substituting the approximations in Equations 93 - 98 into Equations 71
and 72 and summing over the elements gives the discrete system of equations:
44
- A‘l<U^iUj #i>n., + <V5,U5 *,>„J
- Atl<(jL)kD„„k + <(jL)kD„,k 01>()j]
i=l, ...N
J.. % ^ <f'jUi*Vi.
= <(1 - t^rij) Vj0j.
- ^jEhiAt [<VkH|i |^^)>(i^,+ <VjH|, ^^)>nj
- + Ci“ - “^“) + (?!i + <i - ^)1 *i>».,
+ l^t(^)j“ + (^)il *l>nel
- At(<U5iVj[ 4k,>„^,+ <V‘,Vk ^j>„J
- At[<(g^)'‘D„vJ «i>n^j + <(h^)'‘Dvv5 ^i>nj]
i=l, ...N
Equations 99 and 100 can be rewritten as:
“ ®el r At . . m . . At . I
_E__ E [(1 + ^4j)M(!)uf •■ - + ^,Et,AtMS?)u5'
= (1 - firjJjM^pUV + Aif.kjjWvk - ,3,Et,AtM!3)u5
- + (5“ - a^-) + (^ + C$ - o»^)l
+ + vJimWuJ)
- At(^)'‘lM(pD„„f + mWDuv}]] i=l. ...N
(100)
(101)
I?. 7^^il)MipV5- + fif'jMSjJuJ*' + ^.E,,AtM(fvJ-
45
58. The M(p matrices on the left side of Equations 101 and 102 and in the
first two terms on the right side of these equations are lumped so that all elements
are added onto the diagonal. The matrices on the left side of Equations 101
and 102 are decomposed into diagonal and non-diagonal matrices. The non-diagonal
portion of is moved to the right side of the equations. These operations give:
- + <j“ -
+ - At(US,M(I)u{ + Vj,MSf)up
- At(Hl)'‘lMWD„„} + MWD„,y]] i=l, ...N (103)
and
E [(1 +
C 1 *1 J ” ^ ^
= (1 - |irt5)Ml5‘')vy - - E^5At(^,M(f'>)vV->+ ^jMS^V})
- + (}•' - + (^ + Cf - ao})]
+ + (^)jl - + v5,m(«)v})
- Al(3L)k[M5pD„vf + Mj?)Dvv5l] i=l, ...N (104)
where
= the diagonally lumped elemental matrix
- the diagonal portion of the elemental matrix
46
Mlf*) = the non-diagonal portion of the elemental matrix
59. The fully discretized momentum equations can be written in compact form
M Uel
E S [Mtn - [Mfi ] {VVn = E {Pt }
1=1 j =1 J ‘ J ■* J J -> el=l
i=l, ...N
M Uel
E E [Mfi 1 ] {VV‘} = 2 {P; } i=l, ...N
1 *1 j =1 •' ■* ■' •* e 1 *1
where
Mjf = (1 + + /3.Et,AtM(f )
j^2me
ij 2 J
pI»E ^ "£‘ ^
- EuAt(/?,M(f + ^kHWuf)
- 8|iM(T)[(^' + CV“ - m/J-) + (^ + Cf - 0^))
+ + (gt)S - ^‘(''^1 + V‘,M}?)uf)
- At(gL)'‘(M(I)D„„J + m(^)d„,J]] i=l, ...N
pT"'= = ”e‘ ^(1 - |iTi5)M(‘‘'V5 -
- EhjAt(/J,M(f "Vf +
- + c;-‘ - 0^") + (^ + C{ - 0-^)1
+ + (^)j! - luSiwSPvJ + v‘,m(?)v5]
- At(jjL)k[M(pD„v} + m(?)d,vSi] i=l, ...N
j.
+ -r*
i=l, ...N
(105)
(106)
(107a)
(107b)
(107c)
(107d)
In the three-dimensional model options, if the computed friction coefficient exceeds the
maximum allowable value, the bottom stress terms appear explicitly on the right side
XME YMF
of the momentum equations and therefore are included in Pj and P^ :
j7j ^ u
47
- ■ “’^>1
+ + (^)5] - At(Uki M(pU5 + Vk,M(?)uy)
- At(5^)'‘lM(I)D„„5 + m(?)d„v}| -
p^E ^ “|,> - ^m(}‘')u}
- Et,At(/3,M(f”)vf> +
- + <l" - “^*') +(?!+<?-
+ + (^)5| - At(Uj,M(J)v5 + V‘,M(pv5)
- ^‘(h^)‘1m1?D»vJ + M(f)Dv.;] - AtM(?)(^);'
(108a)
(108b)
60. Global assembly and enforcement of the C° functional continuity
requirement gives the following systems of equations:
N
S [eM|f] {gUf ‘} - {gM?f] {gV^*‘} = i=l, ...N
j j j
N
E {gUj*‘} + [gM{^®] {gyf = {gpf**^} i=l, ...N
where
gM^^®, = global diagonal system matrices
gpXME^ gp^E _ giQij^ right-hand-side load vectors,
gUj*^ gVj*‘ = global velocity vectors in the x and y directions
(109)
(110)
Solution StrateE
61. The horizontal discretization for ADCIRC has been implemented with
three-node linear triangles and four-node bilinear quadrilaterals. The triangle element
provides a maximum degree of flexibility and is extremely cost-effective on a per-node
basis for long wave computations. All elemental matrices through are
integrated using a numerical quadrature rule that is specified with the input data. A
48
four-point Gaussian quadrature rule integrates the elemental matrices, Equation 87,
exactly (Connor and Brebbia 1976). However, for most applications, a three-point
Gaussian quadrature rule appears to be sufQcient. The elemental are
computed once and then stored for use during the time-stepping operations.
62. The GWCE is solved first. The global system matrix for the GWCE,
SMjj , is time-independent and is therefore assembled and LU decomposed only
once. ®Mij has a banded structure with a band width that is dependent on the
node numbering of the grid. Prior to running ADCIRC, the node numbering should
be optimized to minimize the maximum difference in the global node numbers
GWCE
associated with each element in the grid. The right side load vector, sp. ^ ig
efficiently updated every time step since all elemental matrices have been
pre-computed. Flux-spedfied boundary conditions have been incorporated into the
load vector by the model formulation and therefore require no additional equation
manipulation. Elevation-specified boundary conditions are incorporated into the
system matrix, , by zeroing out rows corresponding to boundary nodes with
specified elevations and placing a value of unity onto the diagonal. The boundary
condition values are then stored into the appropriate location in The
equations associated with elevation-specified boundary conditions are multiplied through
by a constant to ensure that the modified global matrix has a good condition number.
The modified system of equations (i.e., Equation 92 modified to include the elevation-
specified boundary conditions), is then solved for elevation at the new time level, k+1.
63. The momentum equations are solved second and use the elevation values
at time level k+1 computed from the GWCE. The global system matrices for the
momentum equations, and ), are time-dependent and therefore need to
be reevaluated at every time step. However, since these matrices are diagonal, matrix
evaluation and decomposition are very economical. Specified normal flux boundary
conditions are incorporated into Equations 109 and 110 by reorienting the x and y
equation pairs that correspond to the specified flux boundary nodes into a locally (for
each node) normal/tangential coordinate system. The i3oriented equations are then
replaced by the corresponding specified normal flow boundary condition values (Wang
and Connor 1975; Gray 1984).
64. The right sides of Equations 109 and 110 are dependent on Uj,
and Vj, which are all known quantities, and on and (because of the lateral
49
closure model). Therefore, these equations must be solved iteratively for velodti^ at
the new time level, k+1. and sp. are updated each iteration with the new
values of Uj** and Vj*‘ until a specified convergence criteria has been reached. When
is zero, no iteration takes place.
65. When a three-dimensional option is used, the external mode solution
depends on the internal mode solution through Duu) Duv> Dw, Cf and 7 [or if the
drag coefficients exceed the maximum allowable values on and These
quantities are computed at each internal mode time step and assumed to be constant
in time for subsequent external mode time steps. If the external mode solution and
the internal mode solution are evaluated at the same model time, the external mode
solution is evaluated first. The updated surface elevations and depth-averaged
velocities are then used in the internal mode solution. This solution sequence requires
the specification of initial values i c Duu> Duv, Dw, Cf, and 7 as input parameters for
the external mode solution.
Fourier Properties of the External Mode Solution
66. Fourier analysis characterizes the damping and phase propagation properties
of a numerical solution in relation to the corresponding analytical solution. Although
it is typically applied to the one-dimensional form of the shallow-water equations and
a constant bathymetric depth is usually assumed, the results give a good indication of
how a circulation model will behave in a more general two-dimensional, nonlinear field
application. They also allow inter-comparisons with other discretization strategies.
Procedures for applying Fourier analysis to the shallow-water equations are described
by Pinder and Gray (1977) and Lynch (1978).
67. The discrete form of the ADCIRC 2DDI governing equations has been
Fourier analyzed. These results are presented below along with results from the
Fourier analyses of several other numerical solution schemes for the shallow-water
equations. All of the other numerical schemes that were considered use primitive
formulations of the shallow-water equations (as opposed to the generalized wave-
continuity formulation used in ADCIRC). The schemes include a finite element
solution using linear elements (PEFE) (Wang and Connor 1975; Westerink, Connor
and Stolzenbach 1987, 1988), a second order, non-staggered, finite difference solution
(PENSFD), and a second order, staggered, finite difference solution (PESFD) (Hansen
50
1956; Leendertse 1967). A second order, Crank-Nicholson scheme was used to
integrate the PEFE, PENSFD, and PESFD in time. As described previously,
ADCIRC uses a three-time-level scheme for the GWCE (oi = 0.35, oj = 0.30 and
as = 0.35), and a Crank-Nicholson scheme for the momentum equations. The bottom
friction coefficient (Equation 46) in each model was specified as
r* = 0.8 If /g5~/A (111)
where X is the wavelength of the Fourier component.
68. The modulus of the propagation factor indicates the ratio of the numerical
amplitude to the analytical amplitude during the propagation of one wavelength. The
phase of the propagation factor indicates the phase lag or lead a given wavelength
experiences during one period. Figure 1 pr^ents the modulus and phase of the
propagation factor for the PEFE and PENSFD schemes. Comparisons are shown for
Cr = 0.1, 0.5, 1.0, and 2.0 where Cr is the Courant number based on wave celerity.
At is the time step, and Ax is the gnd spacing. For increasing Cr, both the PEFE
and the PENSFD solutions have less damping than the analytical solution for low
ratios of A/ Ax. Neither solution, regardless of Cr, propagates energy at the shortest
resolvable wavelength, A = 2Ax. This characteristic of PEFE and PENSFD solutions
accounts for the severe 2Ax numerical noise problems encountered using these schemes.
69. Figure 2 presents the modulus and phase of the propagation factor for the
PESFD scheme and the generalized wave-continuity equation finite element
(GWCEFE) scheme used in ALCIRC. For low ratios of A/Ax, both schemes provide
less damping than the analytical solution and show poorer phase propagation behavior
as Cr increases. For a fixed Cr and A/ Ax, the PESFD scheme has slightly better
damping characteristics, while the GWCEFE scheme has better phase propagation
characteristics. At low Cr, the GWCEFE solution leads the analytical solution. As Cr
increases, the GWCEFE phase propagation factor swings through a zero value
(corresponding to perfect phase behavior) and then develops a phase lag. This
indicates that there will be a local minimum in the time convergence curve with
optimal accuracy being achieved at Cr » 0.5.
70. The primary difference between numerical solutions using PEFE and
PENSFD schemes and numerical solutions Uoing GWCEFE and PESFD schemes is
that the latter schemes propagate energy at A = 2Ax. Propagation of 2Ax waves
corresponds to a non-folded dispersion relationship and prevents two responses from
51
MODULUS OF PROPAGATION FACTOR MODULUS OF PROPAGATION FACTOR
a. Cr = 0.1
J/Ax
b. Cr = 0.5
Figure 1. Modulus and phase of the propagation factor for PEFE and PENSFD
solutions (Continued)
52
MODULUS OF PROPAGATION FACTOR MODULUS OF PROPAGATION FACTOR
VAX VAX
a. Cr = 0.1
^Ax VAX
b. Cr = 0.5
Figure 2. Modulus and phase of the propagation factor for GWCEFE and PESFD
solutions (Continued)
54
MODULUS OF PROPAGATION FACTOR MODULUS OF PROPAGATION FACTOR
developing to a single forcing frequency, i.e., one physical response at the forcing
wavelength and one numerical response at a wavelength near 2Ax (Platzman 1981).
As a consequence, GWCEFE and PESFD schemes do not have the severe 2Ax noise
problem of the PEFE and PENSFD schemes.
Convergence Properties of the External Mode Solution
71. In order to verify the accuracy of the external mode solution of ADCIRC
and to establish convergence properties in space and time, ADCIRC-2DDI was applied
to a modified form of the quarter annular test problems originally developed and
appMed by Lynch and Gray (1978, 1979) and Gray and Lynch (1979). These two-
dimensional, variable-depth test problems were developed to give insight into a
numerical scheme’s 2Ax oscillations and its ability to propagate longer physical waves.
The original geometry and bathymetry of Lynch and Gray (1978, 1979) were modified
as follows. The arc of the annulus was increased to 135 deg*; the inner radius was
decreased to 125,000 ft; the outer radius was increased to 650,000 ft. The resulting
geometry, with three land boundaries and one open ocean boundary, is shown in
Figure 3. A linearly varying bathymetry was used that increased from 50 ft at the
inner radius to 260 ft at the outer radius and a quadratically varying bathymetry was
used that increased from 50 ft at the inner radius to 1,352 ft at the outer radius.
These modifications accomplish two things. First, the modified domains are more
representative of a coastal region that extends to near or beyond the Continental Shelf
break. (In fact, the geometry and bathymetry are idealized approximations to the
New York Bight.) Second, the numerical difficulty of the test problems is increased.
72. A sequence of four discretizations was considered; a 6- by 8-node
discretization (Ar = 105,000 ft), an 11- by 15-node discretization (Ar = 52,500 ft), a
21- by 29-node discretization (Ar = 26,250 ft), and a 41- by 57-node discretization
(Ar = 13,125 ft). These are shown in Figure 4. Grids consisting of linear triangles
and of bilinear quadrilaterals were tested and gave very similar results. Only the
bilinear quadrilateral results are presented here. For each grid, five different time
steps were applied: At = T112/8, At = TM2/I6, At = TM2/32, At = TM2/64, and
At = Tii2/128 where T11I2 is the M2 tidal forcing period equal to 44,712 seconds.
ADCIRC-2DDI was run in its linear mode with an M2 forcing frequency. Therefore,
*A table of factors for converting non-SI units of measurement to SI (metric) units is
presented on page 6.
56
Figure 3
650.000 ft.
Test problem geometry
57
a. 6- by 8-node grid
b. 11- by 15-node grid
Figure 4. Grids used for the test problem (Continued)
58
d. 41- by 57-node grid
Figure 4. (Concluded)
59
the theoretical model response should have included only an Mj wave. The resolution
of the Mj wave provided by the sequence of grids varied between 17 and 312 nodes
per wavelength for the linearly varying bathymetries and between 17 and 703 nodes
per wavelength for the quadratically varying cases. Thus the M3 wave was always
well-resolved. €r varied between 0.13 and 39 for the linearly varying bathymetries
and between 0.13 and 88 for the quadratically varying cases. ADCIRC-2DDI is
unconditionally stable in its linear mode and therefore permits the use of Cr > 1.
73. All cases were forced at the open ocean boundary using ^ = 1.0 8in(u;|^ t)
where 0;^^= 2t/TiI3 is the Ms forcing frequency. All other forcing mechanisms (i.e.,
tidal potential, free surface wind stress and atmospheric pressure gradients) were set to
zero. The Coriolis and advective terms were also neglected. The bottom friction
coefficient was set to r* = 0.0001 and the value of Tq = 0.0001. All total depths
were set equal to the depth to the geoid.
74. The computations were hot-started using the analytical solution for the
specified geometry, bathymetry, and friction coefficient. The computations were then
run for 10 tidal cycles to allow a dynamically steady-state numerical solution to
develop. The elevation and radial velocity solutions at each node were recorded
during the eleventh tidal cycle and were Fourier decomposed. Typical results are
shown in Figure 5 for the sequence of runs using the coarsest grid and the linearly
varying bathymetry. The figures compare the exact analytical solution to the
maximum and minimum ADCIRC-2DDI solution for all nodes at the same radius.
These plots indicate that there are no spurious 2Ax modes in either the radial or
angular directions.
75. Error measures were calculated from comparisons between the harmonically
decomposed numerical solutions and the analytical solutions. These were defined as:
(113a)
(113b)
(113c)
(113d)
60
SINE COMPONENT OF ELEVATION COSINE COMPONENT OF ELEVATION
o o o o o o
<»> (./») X««|.A (O|P0»
61
R Xio^ft) R XlO^ft)
a. At = TM2/128
Figure 5. Elevation and radial velocity components on the 6- by 8- node grid with linear bath3rmetry
(Sheet 1 of 5)
SINE COMPONENT OF ELEVATION COSINE COMPONENT OF ELEVATION
c. At = TM2/32
SINE COMPONENT OF ELEVATION COSINE COMPONENT OF ELEVATION
FigTJie 5. (Sheet 5 of 5)
where
Np = the number of nodes within the grid
= amplitude of the sine component of the analytical elevation solution at node i
bf = amplitude of the cosine component of the analytical elevation solution at node i
n
a^ =: amplitude of the sine component of the numerical elevation solution at node i
amplitude of the sine component of the analytical radial velocity solution at
node i
amplitude of the cosine component of the analytical radial velocity solution at
node i
amplitude of the sine component of the numerical radial velocity solution at
node i
amplitude of the cosine component of the numerical radial velocity solution at
node i
These error measures represent the absolute errors in the sine component of the
elevation solution (El), in the cosine component of the elevation solution (E2), in the
sine component of the radial velocity solution (E3), and in the cosine component of
the radial velocity solution (E4).
76. A summary of the error measures computed for all of the test runs is
presented in Table 2. The error measures are plotted against Cr (the average value
for a given grid) for the linear bathymetry test cases in Figures 6 and 7 (Figure 7 is
a blow-up of the low Cr range in Figure 6), and for the quadratic bathymetry test
cases in Figures 8 and 9 (Figure 9 is a blow-up of the low Cr range in Figure 8).
All errors show good spatial convergence; i.e., the more refined the grid, the lower the
error at any Cr- In time, the errors decrease as Cr decreases, until Cr = 0.9 - 1.75
for the linear bathymetries and Cr = 3.5-7 for the quadratic bathymetries. A well-
defined local error minimum exists for all grids within these Courant ranges for both
the sine and cosine components of the elevation and radial velocity solutions. This
local error minimum occurs because the phase of the propagation factor changes from
a phase lead to a phase lag, passing through a region of almost perfect phase
behavior, near Cr » 0.5 (see Figure 2 and associated discussion). Figures 6-9
suggest that the optimal behavior occurs at somewhat higher values of Cr. These
figures were plotted using the average value of Cr for a given grid. However, the
primary errors are generated in the shallow portions of the domain. If the Cr is
66
r«bl« 2. SunMry of Error Hoaturoo for Alt Convorganco Run*
st ty$ m m 9^ m w% m m ^ wy 9*s ^ ^ m »*y
Ul <41 <41 <41 <41
CM 0 0 «* 0
0 ^ 0 0 W
N^ 0 <X CM
PM 0 4SP 0 »
fo r\4 ** VI fx
KKikss
fX K PS4 U1 <X
0 0 «• XI P'1
iiiiisKiisK
.# a- lA CK P14
a-I PA aaT X> PA
VI <A 0 0 0
PX PX PX ^ X4
<0 0 0 0 (A
<0 <a p’w «-
^ A< rx lA O
x> O ^ o o
PX AJ PA •* <0
m tn m m m
<tt <41 144 <M <44
0 0 ^ 0 m
K. 0 rx ^ 0
K< P#^ 4*4 rx V-'
pp^ ppg ^ pa* rx
fSI ^ PX 'X 0
^ ^ 0 INI «-
0 0 0 PA PX
7 0 x> 0 v<
PA Ai 10 <X ^
VI 0 0 PA PX
oJ aaT p.»* |a{ ^
OOOOO
X4 Ul <44 III W
1A «A lA 0 PA
OOOOO
X> VI PA Ka PA
PX PX PX 0
(A lA lA 0 PA
(44 <14 U4 144 <44
0 »“ PA •• 0
A| 0 <1 A. 0
0 «
lA lA (A O PA
141 144 <44 144 <44
PA K O a~ O
A a- ha. fX —
PX PX O Al «—
m m m m
144 <41 <U <44 <44
0 0 P'1 m N
0 » 9> 0 0
0 k! ATI 0 0
p^ p'l 0 P'1 fX
0 VI 0 0 1/1
«14 ^ 0 0 VI
0 0 P'1 PX Ai
(41 <44 144 <y IX
PX 0 fX VI 0
0 0 PA a- P^
0 A< XI
VI 0 PA PX Al
Ul <41 (44 V ^
Sacsic
PA aA PX ^ lA
0 PA pA PA PA
0 0 PX 0 0
IX PX Ai — VI
0 0 0 »A PA
144 <14 <^ <44 111
9SS:iS
0 0 PA •> P^
0 lA 0 PA PA
<44 <14 <14 Ul *V
$• PA 0 0 0
VI 0 0 0 PX
^ 0 PX ^ 0
A A O PA PA
(41 <44 <14 IM <44
32$ tK
PX *<0 A ^ CO
N1 P'4 fX
144 <44 <44 <14 tl4
^ 0 ui 0 VI
P'1 0 0 m P'1
m 0 P'1 P^ Ai
<4t <44 til <44 <44
px m ^ VI 0
0 0 VI P'1 0
0 0 PA AJ PX
<44 U4 <U <44 <44
VI 5> 0 ^ (X
VI 0 PX 0 ^
VI 0 PA PX «X
(44 (44 <44 <44 <41
3 a- rx PA f1
•- PA 0 0
PA PA PA 0 «X
(44 lU 144 (44 <14
0 0 PA 0 VI
0 4«> PX 0 0
0 0 0 PA PX
(44 144 <44 <U (44
lA A> 0 ^ 0
** 0 0 PX A4
0 (A 0 lA PX
<44 (41 <44 (44 <44
O PA a^ O O
VI O ^ P^ PA
A O O PA PX
<44 <41 <44 <41 (44
fXOA4«-a-
O — PA O PA
^ m r>t ^ ^ t>o^ tfy f>t ^ ^ tu r>s t>4 ^•»^>AiMi^
O •- PX PA Pte
PA >0 PX O O
*- PX fA ^ o
O PX O O Aa
N.
(X PA A. O O
PX O CO K. *
(A
Aa ON
PA 4^ O O •
O O A» a CO
• a a O PA
PX o O <
Aa O K O 8
O «** PV A •
PX O 40 —
>a
IS
O A o o
X) A ^ A •
Aa A O O 8
A ^ A a a
|A O O 8 8
IX
o o o
O o O CA •
A • • • ►^
O* O p- PX o
o* •• Al o »
PX O (A ^
X
O a-a PX A
—* PX A ^
PX «4> ^ O
A •- PX O *
PA K> PA K> O
^A4aO«-
O O O ^ PX
K PA K o
PX (A O ^ A4
O O ^ PX O
PA ^ O ^ O
A O a- (X VI
O ^ PX o* o
o
OP
ui
A Aa A Aa
p- PX A O
O O o ^
o
PX
4^ A ^ O K
A4 A O •• PX
O e a* PX O
A 4a. o Aa O
A O •» PX A
O PX O CO
Aa O 4I O S
O^PXA-
— PX 0 0 ^
s33Ss
o 9 o
S2SfcS
Na O 9 O
ssjsfcs
O 9 O
SSSfcS
<1
flK:
3
A
A
g3SSa
^ V DO
SSalfcS
5:838 a
0 •• A A. A
0*000*^
o o o o
0 0*00-^
a o* o o ^
o o o o
o o o o *•
o o o <o
0000^
>0 'O 'O 'C 'O
_ « «9 CO «g «o ' < • • •
00000 0>9k9»0k» . ^ ^ ^ ^ ^
• «.•« t/cmintntn •—*— ^
Nnvc»i»»i»>Rn
k»k.K»n.n. ^ ^ ^ ^ ^ f , , , t
• I . . #i,,f 'O'O'O'O'O
rj fyt fs4 ^ f\t f0y 9n tn tn tn • • • • »
^ ^ ^ ^
^ xT 't# -^r
COflOOCOO S.^«.K.K.»«ia .
O0a^0*0» . .....
^w^J^lJN!KI CCICCIC 2S22S
OCOOOO ^ t • » » •
»«•*« «*««* ••«!> 'O'O'O'O'O
— ^ — fM^M^^»^J^^^ Ro^pii^g^ga^^ .
«*;*3(* SSSSS 332??
99990 00000 V^tf^lAcOm
00000 v%<ninv>y> ^(Viv<voj
m«A«ntnRn r>i f>4 r>i rn
*9999 99999 OOOOO t/^mtfxnin
*9999 OOOOO iOR/«RntnRn ^J<^crM(^lrv
»0000 (OmiAVtlA a-a-a-—*-
«wfvf«4«vfv OOOOO in tn ffy
>0000 Rninic^inm f^#R|f*j«Mf'J
O O O *# O <040000 OOOOO
^ rz ^ in«AM^<n«0
l»*.^a^.^>a^a. JCWJJJIJIC l(JCKXI< IcWvW^
H H K K M OOOOO OOOOO OOOOO
inv^mmi/n fMfM<^r<A«v
OOOOO OOOOO
nifyoirvcv
KKMXIC MJ«KieK
XXXXX OOOOO OOOOO
«nmt/>mi/% (MrMCMrviRV
00000
tn m tn m m
X X X X X
OOOOO
OOOOO
67
SINE COMPONENT OF ELEVATION COSINE COMPONENT OF ELEVATION
Figure 6. Convergence properties of the linear bathymetry test cases
SINE COMPONENT OF ELEVATION COSINE COMPONENT OF ELEVATION
Figure 7. Detail of the convergence properties of the linear bathymetry test cases at low C;
SINE COMPONENT OF EL? ATIO^' COSINE COMPONENT OF ELEVATION
Figure 8. Convergence properties of the quadratic bathymetry test
SINE COMPONENT OF ELEVATION COSINE COMPONENT OF ELEVATION
convergence properties of the quadratic bathymetry test cases at low C
adjusted to account for this, the optimal range of values changes to Cr = 0.52 - 1.07
for the linearly varying bathymetries and Cr = 1-07 - 2.14 for the quadratically
varying bathymetries.
77. It is concluded that the external mode solution used in ADCIRC has
excellent numerical properties. There are no spurious 2Ax or 2At modes due to the
ability of the GWCEFE scheme to propagate high wave number energy. Convergence
properties in space and time are good with superconvergence occurring in the range
Cr = 0.5 - 1.5. In this range, more accurate solutions are obtained using larger time
steps.
Application of ADCIRC-2DDI to the English Channel and Southern North Sea
78. The accuracy and behavio’^al characteristics of the external mode solution
have been tested in field applications including (a) tidal and hurricane storm surge
simulations in the Gulf of Mexico (Westerink et al., in review), (b) tidal simulations
in the English Channel and Southern North Sea, (c) tidal simulations in a small
coastal inlet (Luettich, Birkhahn, and Westerink 1991) and (d) tidal simulations in
the New York Bight. The English Channel/Southern North Sea system is probably
the best documented field site presently in existence for testing a long-wave,
hydrodynamic model. Since the emphasis of this report is on the development and
testing of the various components of ADCIRC, the results of applying ADCIRC-2DDI
to the English Channel and Southern North Sea are presented below.
79. In the mid-1980’s considerable effort was put forth to establish and make
readily available a set of standard grids, boundary conditions, and verification data for
model evaluation in the English Channel and Southern North Sea (Werner and Lynch
1988). This data has been used as the basis for modeling studies for the Tidal Flow
Forum I at the Conference on Finite Elements in Water Resources, Lisbon, Portugal,
in 1986 and for the Tidal Flow Forum II at the VII International Conference on
Computational Methods in Water Resources, Cambridge, MA in 1988. Two collections
of scientific papers have been published from this work and can be found in Advances
in Water Resources. Vol. 10, No. 3 (1987) and Advances in Water Resources. Vol. 12,
Nos. 3 and 4 (Dec 1989).
80. The fully nonlinear version of ADCIRC-2DDI was applied to the grid and
bathymetry shown in Figure 10. The grid consists of 990 nodes and 1,762 linear
triangular elements. The model was forced by specifying 11 harmonic constituents for
elevation (Oi, K|, M2, N2, S2, K2, MS4, MN4, M4, Me, 2MSa) along the two open
72
73
modd boundaries. Wind stress and tidal potential forcings were not used in the
modd runs. Modd parameters were selected to match those used by previous
investigators to allow the direct comparison of model results with field data and with
previously published modd results. The following parameter values were used in the
model: u = 0.0002s-S At = 360s, Cf = 0.002322, f = 0.0001 13341s ‘i, and Ehj = 0.0.
The time integration coeffidents in the GWCE were set to Oi = 0.35, 02 = 0.30, and
03 = 0.35.
81. ADCIRC-2DDI was run for the short-term test case suggested by Werner
and Lynch (1988) covering the period from 0 hr on 15 March 1976 to 24 hr on 17
March 1976. Werner and Lynch (1987) found that it was necessary to use a
minimum bathymetric depth of 15 m throughout the model domain to avoid
generating negative water depths during their simulations. ADCIRC-2DDI ran
successfully using a minimum bathymetric depth as smaU as 10 m, although the
simulated results were highly insensitive to this change at the 19 locations where
observational data were available (see Figures 11 and 12).
82. The first 47 hr 10 min of the simulation were used as a transient start-up
period. Figures 11 and 12 present comparisons between modeled time series and
observed time series of free surface elevation (at 11 stations) and depth-averaged
current speed and direction (at 8 stations) for the final 24 hr 50 min of the
simulation. (The locations of the elevation and velocity stations are shown in
Figure 10. The observed time series were actually reconstructed from 11 primary tidal
constituents at each station. The tidal constituents correspond to those used to force
the model open boundaries and were extracted from raw time series at each
observation station using harmonic analysis.) In general, the model does a good job
of simulating the observed results. Some of the differences can be attributed to local
topographic and bathymetric effects and to the inherent problems associated with
representing bottom stress in a depth-integrated model. Also, Werner and Lynch
(1989) point out that the model results contain harmonic constituents, generated by
nonlinear interactions within the domain, that are not included in the reconstructed
observed time series. By filtering this energy out of the model results, they were able
to reduce the average difference between the simulated and observed surface elevations
by approximately 40 percent. The worst comparison occurs at the tidal elevation
station at Christchurch and is at least partially due to the neglect of the channel
between the Isle of Wight and the mainland (located approximately 25 km east of
Christchurch) in the model grid.
83. ADCIRC-2DDI was also run for the long-term test case suggested by
74
observation - - 15m minimum depth . 10m minimum depth
Figure 11. (Sheet 2 of 3)
observation - - 15m minimum depth . 10m minimum depth
Figure 11. (Sheet 3 of 3)
15m minimum depth . 10m minimum depth
(oas/uj) paada uoipauip
(aas/uj) paads uoipajjp
cx
ao
CN
«o
o
o
o*
ao
CN
(O
o
78
time (hre)
Figure 12. Time series of simulated and observed tidal speed aud direction at various stations (Continued)
Werner and Lynch (1988) covering the 190-day period starting at 0 hr on 15 March
1976. The first 5 days were discarded to allow for start-up transients and the
remaining 185 days were harmonically analyzed using the least squares package of
Foreman (1977). The amplitudes and phases of the primary surface elevation
constituents from the ADCIRC simulation, from a simulation by Werner and Lynch
(1989), and from the observed time series at the 11 elevation stations are compared in
Table 3. The overall comparison between model results and observations is reasonable
considering no effort has been made to calibrate the model by adjusting the bottom
friction coefficient, as attempted by Baptista, Westerink, and Turner (1989). Some of
the largest differences in phase occur at stations that are close to amphidromes. This
is because a small displacement of an amphidrome’s position can result in a large
change in the nearby phase values. Some of the largest relative differences in
amplitude (i.e., percent difference between the simulated and observed amplitude)
occur in the Ma constituents. Bottom friction is the primary nonlinear generating
mechanism for this constituent, suggesting that this process is not captured very well
by a depth-integrated model.
84. Figure 13 presents co-tidal charts for the entire domain for 14 tidal
constituents. The ADCIRC-2DDI results presented in Figure 13 and Table 3 compare
very closely with those of Werner and Lynch (1989). This is expected since Werner
and Lynch (1989) used a depth-integrated, finite element, GWCE-based model that is
similar to ADCIRC-2DDI. The minor deviations between the models are due to
ADCIRC’s use of a non-conservative formulation of the advective terms in the GWCE
as well as slight differences in the discretizations of several of the terms. The close
correspondence between the model results provides an excellent verification of the
formulation and numerical discretization used in the external mode of ADCIRC.
80
Titolt 3. MpHtudt (<n M«t»n) tnd Ptiw I>»t»t1v» to Cnwmiidi (lleu<d«d to th» m— r««t Omgrtr) for th» fit>d D«t«, lwa«l
« «j5§
\ S§i
Mil m
ILl m IT) tn
ft Kt Kt fir)
* ssa
V «ss
« Cti
fi
s
<
o o o
s ass
' 11?
ggg
n§g §^s
SCR S'*-'
fW fNT
ftSI
s?? n»n
o a o o o o
j("” sss
???
8P8
ooo
W.
OOO
g§8
ooo
HP
ooo
eias
|88
)Q S M
•'I
rsi 8<^ Q
ooo
ooo
ooo
HP.
OOO
r*
Hi
SS&
saa
ssr
N 72 fs«
ooo
C P M
OOO
Ui
ooo
§55
ooo
aaa
O' O tA
7 *«
ni fw fs»
0 6/9 80
^ 80
Hi
aAi<
M «
ooo
O ->4 t/%
tqpin
ooo
ooo
HP
ooo
Hi
Mi
t>tO- to
p g a
152
ooo
ooo
P.a
ooo
P.l
ooo
M *3 8^
«0 Kt Kl
M Kl
diCR:
K O
ooo
OOO
OOO
HP
ooo
S ?3 rl
S S S
eias
N fV
ooo
P.s
ooo
safe
'»•»*#
ooo*
P -» K»
£RjR|
ooo
55S
SSS
«♦ 689 6rt
^ m
1/9 ¥> 6/6
^4 fsi
553
rM ra (M
o rst 8<'»
asa
PH
OOO
O 9-^ 9-
#'•1
a5s
tn
«n o> o
nf^9^
Rxia
S9?
ooo
i K>80
OOO
D r'> o
saa
ooo
ooo
*- 9- rsi
ooo
^ O'* o
9^m m
ass
o o
•9 1/9 «/6
O
AAA
ooo
ooo
188
ooo
EER
OOO
ooo
l::8S
K O o
^ «
5?l
SPR
r« K. « o> -4
rM o o 7
§11 i§l
o o o o o o
SDix di^
s^s esa
o S o ®
o o o o o o
a^a 538
8?" ?22
8?? ?22 §2S
ooo ooo ooo
2^1 gas
ooo ooo
««a a-”
§sa 8gg
ooo ooo
aaa a""
cSi sas
9?? saa 355
ooo ooo ooo
91? 855
128 5??
??§ §S2 128 55?
o’oo 0*00 ooo ooo
B S
r|K ||K -;|l.» Ijis llh ip> .{is el!! Il
Isi llJi Ilii iis5 iis3 Ssss shs iiss il-'S si
Q.140 1M. 0.980 521. 0.150 234. 0.680 257. 0.200 296. 0.050 296. 0.020 311. 0.040 334 . 0.040 23 . 0.040 116 . 0.040 160.
0.135 160 . 0.114 334 . 0.116 232. 0.602 266. 0.195 502. 0.052 304. 0.015 324 . 0.038 318. 0.044 21. 0.026 179 . 0.049 257.
0.136 160 . 0.115 354 . 0.116 233 . 0.610 266 . 0.196 302. 0.052 305 . 0.016 329. 0.037 325. O.OU 23 . 0.028 174 . 0.053 252.
a. Oi constituent
Figure 13. Co-tidai charts for simnlated constituents (Sheet 1 of 7)
82
c. Mj constituent
d. Nj constituent
Figure 13. (Sheet 2 of 7)
83
e. Sj constitueiit
f. Ka constituent
Figure 13. (Sheet 3 of 7)
g. 2MS2 constituent
h. 2MN2 constituent
Figure 13. (Sheet 4 of 7)
Ic. MS4 constitneiit
1. Me constituent
Figure 13. (Sheet 6 of 7)
m. 2MSe constituent
n. Residual constituent
Figure 13. (Sheet 7 of 7)
88
PART IV: INTERNAL MODE SOLUTION
Definition and ApplicabilitY of a 3DL Model
85. As discussed in Part II, mode splitting replaces the direct solution of the
three-dimensional governing equations with an ”external mode” computation for free
surface elevation (using the vertically integrated governing equations) and an "internal
mode" computation for the vertical profile of velocity. It was noted in Part II that
all of the physics contained in the three-dimensional governing equations are included
in the vertically integrated equations if the bottom stress and the momentum
dispersion terms are specified correctly. The simple parameterizations of bottom stress
and momentum diSj,>ersion in terms of depth-averaged velocity (Equations 28 - 31) are
physically correct only for the simplest flows (e.g., a logarithmic velocity profile over
depth). Mode splitting replaces these simple parameterizations with the internal mode
equations. Therefore, when the complete internal mode equations are solved, the
bottom stress and momentum dispersion used in the vertically integrated equations are
(in theory) completely consistent with the three-dimensional equations.
86. While the external mode equations are two-dimensional, the internal mode
equations retain the spatial variation of velocity in three dimensions. Considerable
computational savings can be realized if the advective terms and the horizontal
momentum diffusion terms are dropped in the internal mode computations (Nihoul and
Djenidi 1987; Davies 1988). This simplification eliminates all horizontal gradients from
the internal mode equations, thereby reducing them to one-dimensional equations in
space (over the vertical). When simplified internal mode equations are solved, the
bottom stress and momentum dispersion are no longer completely consistent with the
three-dimensional equations. However, these approximations should be physically
correct for flows in which the vertical distribution of momentum at each horizontal
grid point is determined by a local balance between the surface and bottom stresses,
vertical momentum diffusion, the Coriolis force, and the local inertia. (Clearly, this
should encompass a much wider range of flows than parameterizations solely in terms
of the depth-averaged velocity.) The required balance will exist when the rate of
vertical momentum transport is much greater than the rate of horizontal momentum
transport. Assuming horizontal momentum transport is dominated by advection, the
rate of vertical momentum transport will be much greater than the rate of horizontal
momentum transport in the three-dimensional governing equations if
89
» u
da
m
Scaling this yields
EvcUc
»
Ug
Lc ^ ^ hcUc
or TT" ~T? —
He VC
where Eye, Uc, he, and Lc are a characteristic vertical eddy viscosity, horizontal
velocity, water depth, and horizontal length scale, respectively. Dimensional arguments
suggest Eve » ^cUc where ^ is a constant whose value for tidal and wind-driven
flows typically ranges from 10‘3 to 10*2 (Bowden, Fairmairn, and Huges 1959; Csanady
1976; Fischer et al. 1979; Davies 1985). Therefore, the simplified internal mode
equations should be an accurate approximation to the full internal mode equations
provided
^ » 100 - 1,000
Since coastal and shelf waters are usually characterized by large length-to-depth
scales, a model based on the simplified internal mode equations should be widely
applicable in these waters.
87. The model based on the simplified internal mode equations will be called a
three-dimensional local (SDL) model. This name emphasizes the fact that the
simplified internal mode equations give values of bottom stress and momentum
dispersion for the two-dimensional (external mode) equations that are not fully
consistent with three-dimensional equations, but rather are based on a local
approximation of the three-dimensional equations.
Rationale for the DSS Technique
88- Despite the savings gained by simplifying the internal mode equations in
the SDL model, the internal mode equations are difficult to solve numerically because
of the high velocity gradients that characterize the water column near the bottom and
surface boundaries and across strong density changes. Existing state-of-the-art
circulation models use velocity as a dependent variable and therefore require a fine
numerical discretization to resolve regions of rapid velocity change. Davies (1991) and
Davies and Jones (1991) have examined the computational effort required to resolve a
bottom boundary layer using a one-dimensional model through the vertical solved with
finite differences and several coordinate transformation/grid stretching schemes. For
90
tidal flows having an eddy viscosity that is constant over the upper 80 percent of the
water column and that varies linearly with distance from the bed over the bottom 20
percent, Davies (1991) found that it was necessary to use a logarithmic or log-linear
coordinate transformation and at least 20 grid cells to obtain convergence of the
velocity sclation. When the eddy viscosity was determined using a level 2-1/2
turbulent closure, the most efficient solution was found to require a log^inear
coordinate transformation and 50 - 100 grid cells over the vertical for both a
turbulent kinetic energy transport equation and the momentum equation.
89. Practical geophysical flows often have two or more regions containing sharp
velocity gradients over the vertical. Because of the computational overhead in time
and memory required to resolve these features, existing multi-<limensionai circulation
models almost always omit the near bottom region and use a slip boundary condition
that expresses bottom stress as a quadratic function of near bottom velocity. This
assumption is physically correct only when the velocity profile below the lowest grid
point is logarithmic. An accurate treatment of surface and/or internal boundary
layers requires a fine grid in the regions of these layers. In many cases the required
computational overhead makes it impractical to resolve these features in multi¬
dimensional computations. A survey of the recent literature suggests that only rarely
have more than ~ 20 grid cells been used over the vertical in three-dimensional
engineering or geophysical model applications. For example, Oey, Mellor, and Hires
(1985^ used 11 grid cells over the vertical in their model of the Hudson-Raritan
Estuary. Clearly, such models have limited ability to resolve even one significantly
sheared velocity gradient region. (Note: Davies and Jones (1990) have recently
published results from a three-dimensional model of the northern European continental
shelf using 45 grid cells over the depth. However, this model uses a coarse horizontal
grid and omits the advective terms in both the internal and the external mode-
governing equations.)
90. It is well-established from laboratory and field experiments, theoretical
arguments, and conventional one-dimensional models that the time-averaged vertical
shear stress varies rather smoothly through the water column, particularly near
boundaries. Therefore, it should be possible to use a relatively coarse vertical
discretization to solve numerically for the vertical shear stress, even in boundary
layers. A novel technique has been developed that allows the vertical shear stress to
be used in place of velocity as the dependent variable in the internal mode equations.
Applications of the DSS technique using linearized equations of motion (discussed in
detail below) have shown that it provides a highly efficient means of solving the
91
internal mode equations. This technique promises to be invaluable for modeling
coastal and shdf circulation in which the bottom and surface boundary layers comprise
a significant portion of the water column and for modeling processes that are critically
dependent on boundary layer physics such as wave-current interaction, sediment
transport, oil spill movement, ice floe movement, energy dissipation, physical-biological
couplings, etc.
Devdopment and Testing of DSS Method No. 1
91. Internal mode equations can be generated by subtracting the vertically
integrated equations firom the three-dimensional equations (Wang 1982; Sheng 1983;
Davies 1985). Using the three-dimensional equations in the a coordinate system
(Equations 19 - 21), the non-conservative vertically integrated momentum equations
(Equations 25 and 26), assuming a constant density fluid, and neglecting advection
and horizontal momentum diffusion terms, the resulting internal mode equations are
92. Using the eddy viscosity rdationships (Equation 34) to express Tzx and r^y
in terms of velocity and either the slip or the no-slip boundary condition (Equation
10) at the bottom. Equations 114 and 115 can be cast entirely in terms of velocity.
Numerical solutions can then be sought for the dependent variables, u and v. This is
the standard velocity solution (VS) approach.
93. Alternatively, Equation 34 can be inverted to obtain expressions for
velocity in terms of stress
» = - u +
b
V = Vb - V +
b
In Equation 116 the definitions of u and v have been used and nonzero slip velocitira
Ub and Vb have been included for generality. Relating Ub and vb to the bottom
stresses, rb* and Tby, via the slip conditions
92
rbx/Po = k Ub = k(ab + U)
Tby/po = k Vb = k(vb + V)
(117a)
(117b)
Equation 116 can be written as
= Ibx _ u + I da
Po^ (a— bj J PqEv
b
' = ^ I ^
(118a)
(118b)
(For a no-slip boundary condition, the terms and do not appear in Equation
118. The no-slip condition is approached as k-^oo.) Substituting Equation 118 into
Equations 114 and 115 gives:
] =
jry Ts
FT
(119)
d H
5T 1 J ^
b b
“ + ^y] = It ^ ^ ^^20)
Equations 119 and 120 have tzx and Tzy as dependent variables and will be called the
DSS‘ internal mode equations. (The superscript 1 is used to identify DSS method
No. 1.) These equations are forced by the external mode solution (U, V, and H) and
the applied surface stress.
94. Equations 119 and 120 contain both integral and differential terms;
therefore, they are well-suited for a spatial discretization in which Tzx and Tzy are
expressed in terms of assumed shape functions such as the spectral or finite element
methods. Depending on the choice of the shape functions and the functional variation
of Ey over the depth, the velocity profile can be recovered from the stress profile by
solving Equation 118 in closed form. This is an important convenience because it
avoids the troublesome operation of numerically integrating the near^ogarithmic
singularity that occurs in Equation 118 when Ey varies with distance from a
boundary. The restrictions that a closed-form solution for Equation 118 impose on
T’zx) Tzy, and Ey are not severe. For example, Tzx, Tzy, and Ey may be expressed in
terms of polynomials that span the vertical globally or in a piecewise manner.
Polynomial variations of Tzx and r^y are consistent with either the spectral method or
the finite element method; for most practical problems, Ey can be approximated as
piecewise linear over the vertical (Fumes 1983; Davies 1987; Chu, Liou, and Flenniken
1989; Jenter and Madsen 1989).
95. The effectiveness of the DSS> technique is evaluated using a simple test
case consisting of flow generated by a specified surface stress aligned in the x-direction
in a wide, straight channel of constant depth with no Coriolis force. An analytical
solution can be found for the linear version of this problem and provides a benchmark
for the numerical solutions. For convenience, the linear governing equations are
repeated below:
dU
(121)
(122)
(123)
(124)
Equations 121 and 122 are the depth-integrated (external mode) continuity and
momentum equations; Equation 123 is the VS internal mode equation; Equation 124 is
the DSS‘ internal mode equation. Since there is no motion in the y-direction, the
y-direction equations and the subscript "x" in the stress terms have been dropped.
96. The Galerkin-spectral method, with shape functions consisting of Legendre
polynomials (LPs) over the interval -1 < cr < 1 is used to discretize the VS and the
DSS* internal mode equations. The order LP is denoted Lb and can be computed
from the recursion formulas
Lo(£r) = 1
Li{a) = ff
= [^] fir - Lr-i
The first eight LPs are shown in Figure 14. Other properties of LPs of note are
Lrs = Lr(l) = 1
Lrb = Lr('l) = (-I)--
d(r = 2
94
Polynomials
I^Lr(o^) d<r = 0 for ail r > 1
It has been shown for wind-driven circulation that velocity solutions using Legendre
and Chebyshev polynomials yield results of virtually identical accuracy, that these are
highly superior to velocity solutions obtained using expansions of trigonometric
functions, and that these are more accurate than velocity solutions computed with a
second-order finite difference scheme having the same number of degrees of freedom
(Davies and Owen 1979; Davies and Stephens 1983). For further information on the
use of spectral methods in three-dimensional circulation models, the interested reader
is referred to an excellent review by Davies (1987).
97. The Galerkin-spectral discretization for the VS internal mode equation is
obtained by multipipng Equation 123 by the weighting fuiiCtion and integrating
from -1 to 1, i.e.,
4|‘t. u d. - I (l, y.^) d, = - J [a - a] |‘i. d, (125)
-1 -1 -1
Integrating the second term in Equation 125 by parts
Ih
Po
1
Tz
Po ^
(126)
and substituting this into Equation 125 yields;
If j I, u d(r + I j g d<7 = - J j L. der
-i -1 -1
[
-I
T Ts.
Using the definition of the LP, Equation 127 simplifies to
0 = 0
1
If I L„ u d(r + I I ^ |§a da - ^ - La,b
-1 -1
Ik
Po
(127)
m = 0
(128)
m > 1
(129)
Since Lo{a) = 1, the operation that generates Equation 128 is equivalent to
integrating Equation 123 over the depth when m = 0. The identity in Equation 128
occurs because Equation 123, by definition, has no depth-averaged component.
96
98. The final steps in applying the Galerkin-spectial discretization to the VS
are to substitute Equation 34 for Tz/po in Equation 129 (noting that ^ and to
expand u as a series of LPs with time-varying coefficients, ^n(t), i.e.,
R
u(ff,t) = J^(t) Ln((r) (130)
n*t
I 1 I p^] ® ^
n»l -1 n»l -I
Because |'Z/n d(r = 0 for n > 1, the necessary condition |'u da = 0 is identically
satisfied by the spectral solution by using only the n > 1 LP. The solution of
Equation 131 requires a bottom boundary condition (Equation 117). After expanding
u, this becomes
R
I ^nb = - U +
n = t
(132)
99. The Galerkin-spectral discretization for the DSS‘ is obtained by
multiplying Equation 124 by the weighting function Ln{a) and integrating from -1 to
1, i.e.,
lit [f I ^ [t It^^) ■ fe] I da - I j L. |5(^) da -
-1 -1 -1 -1
.It - Efc] 1 '
-1
Lm da
Integrating the stress derivative by parts changes Equation 133 to
[f ^ 1 ‘‘‘'‘‘"J + [e Ie<S)
h d
J ^IJ J
-t -1
Using the definition of the LP, Equation 134 simplifies to
(133)
- E ^ ^ - 1 ^ - [It - ffe] I ^
-1 -1
97
0 = 0
m = 0
(135)
jltUi- + + = 036)
-I -i
The id^tity in Equation 135 occurs because Equation 136 has no depth-averaged
component.
100. The final step in applying the Galerkin-spectral method for the DSS* is
to expand rz/po as a series of LPs with time-varying coefficients, an(t), i.e.,
^ ^ I «»(‘) i«(^)
» 1
» 1
I" I l?^| ^ ^ I J ^nbfims ] -
(137)
m > 1 (138)
n = 0 '-i
n = 0 -I
101. The bottom boundary condition was introduced into Equation 118 and
subsequently into Equation 138. Therefore it does not generate an extra equation, as
was the case for the VS. However, the stress expansion, Equation 137, does not
automatically satisfy the condition that ^ u da = 0. Rather this must be enforced
explicitly. Using Equation 118 and the definition of u, this requirement generates the
additional equation
h
da = 0
(139)
Substituting the expansion for rz/po into Equation 139 yields
(140)
102. The relative merit of the DSS‘ versus the VS was evaluated by
comparing solutions computed numerically with analytical solutions for the problem of
wind-driven circulation in a closed, rectangular channel aligned with the x-axis and
having a constant bathymetric depth. This was done for a steady-state case, for a
periodically varying wind stress, and for an instantaneously imposed wind stress.
103. In each test case, Ey was assumed to be linear over the depth as
expressed by
Ev(a) = Ezo(a+l4-ao) (141)
98
where Cq = 2zo/h is the dimensionless roughness height. It is well known &om
theoretical, laboratory, and field experiments that an eddy viscosity that increases
linearly with distance from a solid boundary realistically reproduces the physics of the
boundary layer near the boundary (Monin and Yaglom 1971; Schlichting 1979; Grant
and Madsen 1986). Despite the fact that this does not hold over the entire depth,
(e.g., it has been suggested that Ey should also increase linearly with distance below
the free surface (Jenter and Madsen 1989)), Equation 141 is used here because it
generates a realistic bottom boundary layer and because it simplifies the analyses of
model results by introducing only two parameters, Ego and Co, into the problem. As
is shown below, the presence of a velocity gradient region at the bottom is sufficient
to illustrate the advantage of the DSS over the VS. In fact, the use of an eddy
viscosity that does not also give a boundary layer at the free surface is a considerable
advantage for the VS, since it eliminates the additional need to reproduce velocity
gradients there.
104. Assuming reasonable ranges for Zq of 0.1 to 10 cm, and for h of 1 to
100 m, suggests values of ao ~ 10** to lO'*. (The combination of Zq = 10 cm and
h = 1 m, which gives Cq ~ 10**, is not considered realistic since Zq is typically
3 to 10 percent of the physical roughness height. In this case the physical roughness
would occupy the entire depth.) Assuming the slope of the variation of Ey with z
* * *
scales with Ub, (Ub = V n)/po ), then Ego ~ Ubh. If Ub varies over the range 0.1 to
10 cm/s, Ego ~ 10*5 to 10* m^/s.
105. Equations 131, 132, 138, and 140 show that the VS and DSS* require the
specification of Ta/po (which is the input forcing) and U. To eliminate the possibility
that errors in the solution for U might affect the comparisons, U was obtained for
each test case from an analytical solution of Equations 121 - 123. As a result, errors
in the VS and DSS* over the vertical do not feed back into the solution for U as
they would if the complete problem was solved numerically.
106. In all of the results presented below, bottom stresses are obtained from
the VS by using computed bottom slip velocities and the linear slip boundary
condition (Equation 117). Comparisons indicated that this method gave more accurate
values of bottom stress than those obtained by evaluating Equation 34 at a = -1.
(A similar conclusion was reached by Gresho, Lee, and Sani (1987).) Velocities are
obtained from the DSS* by solving Equation 118 analytically using the computed stress
profiles.
99
(142)
108. The VS and DSS^ are obtained from Equations 131, 132, 138, and 140 by
dropping the time derivatives, setting U = 0, and considering all other terms to be
constant in time.
109. Figure 15 presents a comparison of vertical profiles of horizontal velocity
for several combinations of K and Co computed from the analytical solution, the DSS*
using 2 LPs and the VS using various numbers of LPs. Equation 143 indicates that
Ok
the analytical solution for stress varies linearly over the depth, regardless of the form
of Ey. This solution can be represented exactly by the DSS‘ using only the n = 0
and n = i LP; therefore the DSS' and the analytical solution in Figure 15 are
identical. Equation 146 indicates that the analytical solution for velocity has a
logarithmic variation over the depth and consequently a potentially sharp gradient
region near the bottom. In Figure 15a the combination of a small K (large amount
of slip) and a large <7© minimizes the gradient region. Over most of the depth the
velocity profile is nearly linear and therefore closely reproduced using a VS with 2 LP.
However, approximately 5 LPs are required to capture the mild velocity gradient near
the bed. In Figure 15b, the same K is used with <7o reduced by two orders of
magnitude. This has the effect of pushing the gradient region closer to the bottom
(i.e., it is equivalent to increasing the depth by a factor of 100 for the same
roughness) and therefore steepening the velocity gradient. Because the velocity profile
100
Solution - DSS 2 LP (o) 1 I — Exoct Scrfution - DSS 2
Vertical profiles of horizontal velocity for the steady state test case for the VS and DSS^
is nearly linear over much of the depth, it is reproduced well by the VS with 2 LPs.
However, near the bed, approximately 10 LPs are required for the V3 to capture the
gradient region. As discussed below, this results in a poor prediction of bottom stress.
110. In Figures 15c and 15d, a high value of K is used, resulting in essentially
no slip at the bottom. For large <To (Figure 15c) a velocity expansion of 10 or more
LPs is required to reproduce this profile. Reducing ao by two orders of magnitude
(Figure 15d) sharpens the profile further, and approximately 20 LPs are required to
capture the velocity profile away from the boundary. Many more are required to
represent the gradient region near the boundary.
111. As noted above, an important reason for using a three-dimensional model
in place of a two-dimensional model is the former’s improved representation of the
bottom stress. However, since stress is proportional to the velocity gradient
(Equation 34), or the bottom velocity (Equation 117), the bottom stress may still be
represented poorly if the gradient region near the bottom is not resolved properly. To
illustrate this problem, a comparison was made between the analytical bottom stress
and computed bottom stresses from the DSS‘ and the VS over the practical range of
K and ao- The DSS* reproduces bottom stress exactly using 2 LPs. On the other
hand, Table 4 presents a summary of the number of LPs requited for the computed
bottom stress using the VS to come within 10 percent of the analytical bottom stress
as a function of K and ao- Clearly, it is computationally practical to use the VS
only for large roughnesses and large amounts of slip, both of which tend to minimize
the velocity gradient at the bottom.
112. Although quite simple, the steady-state case demonstrates the relative
ease with which a DSS can resolve a realistic boundary layer (i.e., no bottom slip and
a linearly varying eddy viscosity) in a hydrodynamic model that explicitly includes the
vertical dimension. In the following examples we evaluate how this highly desirable
capability is affected by unsteady conditions. Only the no-slip case (K = 1,000) is
considered.
113. If a periodic surface stress is assumed of the form rs(t)//7o=
(where u is the forcing frequency and i 5 solutions can be sought to
Equations 121 - 123 that have the form U(t) = Ueiu>*i, u(£r,t) = u(a)ei‘^^ rb(t)/po=
{'Th/ stnd r/(t) = (Note: Tg/po, U, u(«r), rb/po, and are all complex
variables; therefore they may be out of phase with each other.) Substituting these
into Equations 121 - 123 transforms the linear hydrodynamic equations into
102
Table 4
Stftadv-State Bottom Stresses Computed Using Velocity Expansiong
7"b( anal -llx c
K
#LP J
■Hjc anal )
10-2
10-1
3
0.100
10-2
10®
8
0.091
10-2
101
g
0.099
10-2
102
10
0.078
10-2
103
10
0.078
10-3
10-1
8
0.096
10-3
100
21
0.098
10-3
101
24
0.095
10-3
102
24
0.098
10-3
103
24
0.099
10-4
10-1
22
0.100
10-4
lO®
<40
0.192*
10-4
101
<40
0.242*
10-4
102
<40
0.249*
10-4
103
<40
0.249*
10-5
10-1
<40
0.174*
10-5
10®
<40
0.476*
10-5
101
<40
0.602*
10-5
102
<40
0.619*
10-5
103
0.620*
* This is the minimum difference obtained using no more than 40 Legend*o polynomials.
103
(147)
(148)
iu}fj + h = 0
1
•wU = -g (r, - Tb)
iu4 - |- f(»-yv
35
(149)
114. The procedure used to solve Equations 147 - 149 analytically, together
with the linear slip boundary condition, has been presented previously (Lynch and
Officer 1985; Lynch and Werner 1987) and is not repeated here. Rather, the solutions
are given without derivation in Table 5 (Equations 150 - 163).
115. Spectral approximations for the periodic case are generated by expressing
/9n(t) = and an(t) = and substituting these as well as the periodic forms
of u(t), rz(<7,t)/po, Ts(t)/po, nj(t)/po and U(t) into Equations 130 - 132, 137, 138, and
140. The resulting equations for the VS ate
X
u(cr)
n 9 f
ri Ln(
(164)
X
n * 1
.1
X
1
n»l
J
LjuL-,
■i
n d<r
+
n»l
J W
-l
dLju J
W
1
«| o
II
? i^.b
Po j
m > 1 (165)
H
i0,
nsl
f'nb
= -
U +
Th
^0
(166)
and for the DSS‘ are
Po
X
= ^ ®n Ljj
n 9 A
(167)
X
M — V
1
H 1
h^ V .
^ 2
n»0
i-
-1
f in
1 e;
- 1
dffdff +
1
I ®n[f
rt*0 -1
^ d<T
d
+ LnbLmbj
_ Is
Po
m > 1 (168)
X
I
IB A
■‘rt
J ^
- dtrderj
U
- ¥
(169)
n*0 -1-1
116. The periodic solution depends on the dimensionless parameters K and Cq
(as found for the steady-state solution), a dimensionless channel length L' , a
dimensionless frequency fi, and the dimensionless position in the channel x/L. L' is
the ratio of the channel length, L, to the wave length of a shallow-water wave having
period u) (Equation 154). fl is the ratio of the time scale lor momentum to be
104
Table 5
Analytical Solution for the Periodic Test Case*
UqEzo
■Ai(a) - Ai'
, 1
"A 2(a) - ki
hTg/po
B
2(2+(To)
L ^
UEzO _ UqEzO
h TsIpq hrg/pol
1 +
2(2+'^) F-
(l-exp(-AL'))exp(AL'^) - (1-exp (AL'))exp(-AL'^)
UqEzo _ i(lj-72)
bTg/po n
L' = —
exp(AL') - exp{-AL')
vP"
tiizo
A = V i7i - 1
n(Ai+B)
72 =
Ai(-1 )A2 _ A2(-l)'j
B(Ai+B)
B
ao_
2+0-0
Ai(o) = /xi(cr)/i2(l) - /X2(«^)Ai(1)
A2(o^) = /ii(a)^^2(-l) - A2(-l)j - M2(o^)|^/ii(-l) - /fiC-l)!
S - Ml(l)|^y^2(-1) - /*2(-l)j - /t2(l)^/^l(-l) -
fii{a) = ber^[n(a+l+Oo)]‘^2j + i bei^[n(o-+l+Oo)]<^*j
^l2{a) = ker ^[n(o+l+cro)]^2j + i kei^[n(o-+l+o-o)]’/*j
(150)
(151)
(152)
(153)
(154)
(155)
(156)
(157)
(158)
(159)
(160)
(161)
(162)
(163)
*ber, bei, ker, kd are zeroth order Kelvin functions, an overdot ( ) = djda,
an overbar ( ) = j /j da
105
transported through the water depth, h^/Ezo, and the forcing time scale, l/w,
Equation 155. Assuming ranges for u of lO'S sec-i to lO'® sec'i, L of 1 km to
10® km, and h and Ub as given previously, suggests L' ~ 10'® to 10® and fl ~ 10 to
10®. In all cases results are presented for x/L = 0.5, as these are representative of
the behavior throughout the rest of the channel.
117. Figures 16 and 17 present magnitude and phase portraits of the velocity
structure for K = 1,000, L' = 1, and four combinations of ao and n. For the case
fl = 10*1, momentum is transported through the depth in only a fraction of the
forcing period. Figures 16a and 16b and 17a and 17b show that the velocity
magnitude and phase obtained from the DSS^ using 2 LPs are virtually identical to
the analytical solution; therefore, the stress variation is very close to linear over the
depth. This linear stress variation suggests that the momentum balance over the
depth is nearly at steady state and is consistent with the low value of fl. Since
steady state is approached as fl -* 0, the DSS‘ using two LPs gives a highly accurate
solution for fl < 10*i as well. The VS is able to capture the phase change through
the water column with a comparable number of LP to the DSS*. However, as was
the case at steady state, for ao = 10'^, approximately 10 LPs are required to
reproduce the velocity magnitude with an accuracy comparable to the DSS* using
2 LPs. For (To = 10'^, more than 20 LPs are required.
118. For the case fl = 10, the vertical momentum balance is no longer near
steady state; consequently the DSS‘ requires more than 2 LPs to capture the vertical
stress variation. Figures 16c and 16d and 17c and 17d suggest that approximately
4 LPs may be needed by the DSS'. The VS, however, requires at least 10 LPs for
ao = 10**, and more than 20 LPs for ao = 10*^.
119. Figures 18 and 19 compare the amplitude and phase behavior of the
analytical solution for bottom stress with solutions obtained using the DSS* and VS.
These runs were made using a single value of ao = 10'*, but varying fl and L'. The
10^ change in L' has minimal effect in these pictures, indicating that the number of
LPs required for the DSS* or the VS to converge to the analytical solution is only
very weakly dependent on L'. For fl < 1, the DSS* with 2 LPs is nearly identical to
the analytical solution, while for larger fl the number of LPs required by the DSS*
increases to as many as 7 for fl = 10®. Considering the fact that comparable results
using the VS require the use of more than 20 LPs, the DSS’ is computationally quite
superior to the VS for all fl.
120. Although the Coriolis force was omitted from these test cases, the results
can be used to infer whether a DSS will be equally effective when the Coriolis force is
106
Exact
107
I Ts/ph I I T,/ph 1
Pirure 16. Vertical profiles of the horizontal velocity magnitude for the periodic test case for the
^ VS and DSSi
Solution f (0/1 r — Exoct Solution
108
Vertical profiles of the horizontal velocity phase for the periodic test case for the VS and DSS‘
_oi«=.i ssa
- (i^hl)ohvj
((*«»)u - ('*“»)5L)oMV J
Figure 19. Ratio of the spectral to
VS and DSS^
included. The counterparts to Equation 123 for the case in which Coriolis is included
are
iwii - fv = _ r*,) (170)
fii = ((a-b^Ev ^ ^ ^171)
It has been shown (Lynch and Officer 1985) that the linear combinations of u and v
transform Equations 170 and 171 into
i(uH-f)v* - ((a-b^Ev ^ + 1(7^^ _ .r8y)j (172)
i(u>-f)v- - ((a-b^Ey ^ _ 7-^^) _ 1(7.^^ _ 7.gy)j (I73)
121. Equations 172 and 173 show that the vertical structures of v* and v' are
uncoupled and that each is analogous to the structure of u in the absence of the
Coriolis force, except that v* is forced by the frequency w + / and v* is forced by the
frequency w - / Therefore the vertical structures of v* and v* will depend on the
dimensionless frequencies fl* = 0 + <9^ and fl* = Q - respectively, where
3^ = /hVE^. At mid-latitudes, /~ 10'< sec*‘, giving the range of 10'* to 10^.
This yields values for U* and 11' in the same range as 11; consequently the results
shown in Figures 16-19 are also indicative of the performance of the DSS‘ and the VS
when the Coriolis force is included in the governing equations.
122. Analytical solutions can be obtained for the test problem for a transient
forcing by decomposing the forcing into its Fourier components, using the periodic
solutions presented above for each Fourier component and superimposing the resulting
periodic solutions. In this section an illustrative set of results for bottom stress are
presented for the often-used problem of an instantaneously imposed wind on an
initially quiescent channel. Representative values of L = 100 km, h = 50 m,
Oq = 0.01, and Ezo= 0.5 m*/s are used.
123. An instantaneously imposed forcing cannot be represented exactly by a
finite Fourier series; however.
-liiO- = i + 2 y
Tg- steady ^ " x( 2n-l )
(174)
111
gives an approximation to a square wave of period T, as shown in Figure 20. By
selecting T to be larger than the time required for the basin to reach steady state
and considering only the period 0 < t/T < 1, a reasonable representation of an
instantaneously imposed wind can be obtained and used to develop an approximate
analytical solution. Sensitivity analyses indicated that when 50 or more terms were
used in Equation 174, minimal change occurred in the analytical solution of the basin
r^ponse and any change that did occur was limited to times very close to zero (i.e.,
on the order of t/T < 1 percent). Seventy-five terms (N=74) were used in Equation
174 for the solution shown in Figure 20 and the runs presented below.
124. The VS and the DSS' for the transient test case were obtained by
discretizing Equations 123 and 124 in time using a Crank-Nicholson scheme. As
discussed above, the analytical solution for U was used to force these equations,
thereby eliminating any feedback of error from the vertical representation into U.
Figure 21 presents a comparison between bottom stresses obtained analytically and
from the VS and the DSS*. The DSS* with 3 LPs is quite close to the analytical
solution except very near t = 0 (due primarily to the overshoot in the forcing in
Figure 20). Conversely, 15 or more LPs are required for the VS to attain comparable
accuracy. We note that this test case uses (t© at the upper limit of the practical
range and therefore is the easiest case for the VS to capture. For smaller values of
(To, the transient performance of the VS becomes even poorer as suggested by the
steady-state results in Table 4.
125. The results of this section suggest that shear stress can be a highly
efficient substitute for velocity as the dependent variable in the internal mode
equations. For this to be accomplished, it is only necessary that the shear stress and
the vertical gradient of velocity be linked via an eddy viscosity relationship.
Depending on the choice of shape functions and the functional variation of eddy
viscosity over the depth, the velocity profile can be recovered from the stress profile
in closed form. Under these conditions the difficulties associated with numerically
integrating a near-logarithmic singularity are avoided. Most practical problems can be
solved subject to this restriction by allowing a global or piecewise polynomial variation
of Tz and a piecewise linear variation of Eg.
126. One disadvantage with the DSS' is that it yields a fully populated matrix
on the left side of the discretized equations that must be reformed, decomposed, and
solved at every time step if a time-varying eddy viscosity is used. This requires ~
operations to solve for stress and ~ N* operations to extract velocity (using
Equation 118), where N is the number of LPs that are used. Although often only a
112
t (hrs) t (hrs)
^ nmncrical and analytical bottom stress for the transient test case for the
few LPs are required for an accurate solution, as N reaches ~ 10, the computational
attractiveness of the DSS‘ rapidly diminishes in comparison to a VS that only requires
the solution of a banded matrix (e.g., Lynch and Werner 1991). Part of the reason
for the fully populated matrix is due to the spectral method’s use of globally, rather
than locally, defined functions. If Equation 124 is discretized using the finite element
method with linear elements, the left-side matrix is the sum of a triangular plus a
tri-diagonal matrix. This requires ~ operations to solve, where M is the number
of nodes used over the depth. It can be shown that the triangular part of the matrix
arises because of the integral term in Equation 124.
Development and Testing of DSS Method No. 2
127. The solution of a fully populated or near-triangular matrix system can be
avoided by reformulating the DSS internal mode equations to eliminate integral terms
£rom the left side. This can be accomplished by generating internal mode equations
by taking the vertical derivative of the three-dimensional momentum equations rather
than by subtracting the vertically integrated equations from the three-dimensional
equations. The use of internal mode equations derived by taking the vertical
derivative of the three-dimensional equations has been reported by Tee (1979).
Although this report focuses primarily on the simplified internal mode equations for a
constant-density fluid, the derivation of the full internal mode equations is presented
below for completeness.
128. Differentiating the ^-coordinate horizontal momentum equations
(Equations 19 and 20) with respect to a, and substituting Equation 21 for dp/dcr
gives
(Note, this is illustrated for the x-momentum equation only. The y-momentum
equation follows directly.)
Using the eddy viscosity relationship for r^x and rzy (Equation 34) the vertical
gradient of velocity can be expressed in terms of the shear stress as
da _ Etzx
W - EvCa=b)po
115
(176b)
dv _ Hrav
'Ua iSv(a-b)po
Substituting Equation 176 and the expansions
-dBm - ^ W
dfda\ _ d f„du\ 5u dv
^ ■ -g? ^
into Equation 175 gives:
H'^ZX
2v(a— b)po
fHrzv (a-b) d^r^x _ d
iv( a-b)po Hpo dfT^ m
ife] ■
vHTzx
^v( a— b)po
_ ^ r wHtzx 1 . Htzz dv
do j_Ev( a— b)po J Ev( a— b)/Jo
Using the additional expansions
^ r wHTzx 1 _ WH dfTzx\ _
gg[Ev( a-b)poJ (a— b)po
d r Htzx 1 _ H d / T zx\ ,
'gt[Ev(a-b)poJ (a-b)po ^Ev ' ^ I
d \ uHtzx 1 _ uH g ,Tzx\ .
iL- .. _ mi o,rzx\ , Tzx PUn
c5c[Ev( a-bjpoj (a-b)po ^E^^ Ev(a-b)po dx
Hrzv gu gbx , gmx
[a-b)po dy da' da
Tzx \dE , dnU , gvHl
+ gr~ +
ZX
-b)Po di
d vHTzx 1 _ VH g ,Tzx\ .
g5^lE;-(a-b>oJ “ (a=b7^ o5^Er^ +
Tzx dvE
2v(a-b)po W~
Equation 177 can be written in final form as
g Ttzx 1 _ fT'zv (a-b)^ d^Tzx _ _ (a — ^b)fgbx 3mxl
EtJo ~gg5- - cx - - -^J
where Cx represents the contribution of the nonlinear advective terms
= _ ii^FXsilJ _ \T,t^ 1 - f '^zx 1 I Tzx dv T zv gu
e;;^ ^ " e^ ^
Applying the same transformation to the y-momentum equation gives
d T zv I f^zx (a~b)^ d^Tzy _ (a — b)r^v gmv
fa-bir^.
where
,,d \Tzv
.d fr
d \t
(177)
/• I ..V T %y d T xy I T XV ^ T xx dv
'7 - - ®
Introducing the complex shear stress Tz = t^x + irzy (where i = ^ -'1 ),
(178)
(179)
(180)
(181)
116
Equations 178 and 180 can be combined into a single complex equation
c + b + m
d
Tz
(a-b)^ d^Tz
where
C — Cx "I" 3iCy
If + i i"'
m =
ftllx 1 : ftHv
3^ + ‘ 3?
(182)
{183a)
(183b)
(183c)
Because both r® and Ey may vary in time, the discretization of Equation 182 in time
may be facilitated by expanding the leading term as:
d Tz _ 1 d fTz\ _ Tz 1 ^Ey
^[EyPoJ ~ 11^ '^po' Po'E^ M
(184)
Substituting Equation 184 into Equation 182 and multiplying both sides by Ey gives
5Ey
'diSpo' Po^y ^
+ if^ -
Po B ^Pq
- Ev[c + 6 + m]
(185)
129. For the 3DL model, the baroclinic, advective, and horizontal turbulent
momentum terms are assumed to be equal to zero. This leaves
d
'di
4. if
(a-b)2 d^Tz
0
(186a)
(186b)
as simplified DSS^ internal mode equations. (The superscript 2 is used to identify
DSS method No. 2.) We note that for an eddy viscosity that is constant in time,
Equations 186a and 186b have the form of complex diffusion equations for stress.
This provides a physical interpretation for the internal mode equation; i.e, it describes
the turbulent diffusion of stress through the water column.
130. Because of the second derivative term in stress in Equations 186a and
186b, two boundary conditions are needed to solve either equation over the vertical.
The free surface boundary condition is
Tz/po = Ts/po at <r = a
(187)
where Tg is the specified surface stress. A second boundary condition can be
117
generated by requiring the depth average of the internal mode velocity to match the
external mode depth-averaged velocity. From Equation 118 this becomes
+ (5^. ITSe; drf. = u + iv (188)
b b
131. To avoid the fully populated matrices generated by the Galerkin-spectral
method, the DSS® uses the Galerkin-finite element method to discretize the internal
mode equation over the vertical, r^/po is expanded over M-1 depth intervals using
depth-dependent, locally defined basis functions Fr(<r) and complex coefficients 7r(x,y,t)
u
Po “
The Galerkin-finite element forms of Equations 186a and 186b are obtained by
substituting Equation 189 for Tz/po, multiplying each equation by Fr(<T) and
integrating with respect to a over the interval from a to b:
^ f a a * .a
F dn
m*l L
and
u
= 0
n = 1, ...M
I I? + - 7i» I ^ da - 7„
■ = iL i i ''
F F dn
n = 1, ...M
(190a)
= 0
(190b)
132. Linear chapeau functions will be used for Fr(a). The tendency observed
in the DSS* results for stress to become linear over the depth for fi < 1 suggests that
these functions should give a good representation of stress if the element size is
selected so that n® <v 1. (Qe is identical to Q except it is defined using the element
size rather than the total depth.) However, Equations 190a and 190b require a
interpolating basis. To lower this requirement to C°, we integrate by parts:
I
I
Fn
da =
i
-I
dFn dFn
W W
da
(191a)
118
(191b)
j E.fn <i» = Ev(a)f„(a)5^ - E.(b)F„(b)^^ - | |5(EvF.)
for Equations 190a and 190b, respectively.
133. Using linear basis functions, when n = 1 and n = M, the first two terms
in Equations 191a and 191b exactly cancel the integral terms in these equations
making the total diffusion terms equal to zero. However, when 2 < n < M-1, the
first two terms in Equations 191a and 191b are identically zero. Therefore, for
2 < n < M-1, Equations 191a and 191b can be substituted into Equations 190a and
190b to give physically meaningful equations:
V \d\^ rVaFn . - (a-b)2r^dFn
Z ■^t [^“ J 'ETT ^ J J ^ \ 'ST '5a~ ^
^ b b b ''
n = 2, ...M-1 (192a)
" [■ * * 2 a
I [^t^ ^ - 7ai^| ^oi-Pnln(Ev) da + 7B^^ga^ da = 0
n = 2, ...M-1 (192b)
134. The boundary conditions are used to supplement Equations 192a and 192b
when n = 1 or n = M. Equation 187 is used in place of the n = M equation:
Re{7y} = ’■sx/Po and lB{7y} = rgy/po
In place of the n = 1 equation, Equation 188 gives
^f2iE|Ibl + I da dal = U + iV
(193)
(194)
135. Velocity is recovered from stress by solving the discretized version of
Equation 118
» + iv = I 2!!E|M + I I®
(195)
136. Equations 192a or 192b and 193 form a tri-diagonal system; Equation 194
adds a fully populated bottom row to this system. However, only a few extra
computations are required to reduce the system to tri-diagonal. Therefore, the
number of operations required to obtain a solution for stress at each time step scales
with M. Since T^fpo is piece wise linear over the depth, the integrals in Equation
195 can be evaluated analytically for many functional forms of Ey. For most practical
model applications, it can be assumed that Ey has a piece wise linear variation with
depth (Fumes 1983; Chu, Liou, and Flenniken 1989; Jenter and Madsen 1989) This
is physically correct near boundaries and makes the analytical solution of the stress
integrals particularly simple. Using this functional form for Ey, the number of
operations required to analytically extract velocity from stress also scales with M.
137. An initial evaluation of the DSS^ has been made using the same test
problems solved for the DSSi and the VS. For these tests f = 0, Ty= 0, and Ey is
constant in time. Therefore, Equations 192 - 195 are simplified to:
1 ^
= 0 n = 2, ...M-1
r I? I j d(7
= 0 n = 2, ...M-1
7y — TsxIPo
n = M
I 7m
Fjb) . h f f
?ok (a-b)2j J E7
b b
da d(T
= U
n = 1
u = Jtb
n-l
Pak ^ (a=BJ J e::
da
(196a)
(196b)
(197)
(198)
(199)
138. To distinguish between the two internal mode equations, results are
designated as DSS| or DSS| depending on whether they are based on Equation 196a
or Equation 196b, respectively. In all of the results, a specified number of equal-sized
elements was used over the vertical. It may be possible to improve the efficiency of
the DSS2 further using elements that are not equally sized. However, this option has
not yet been investigated completely.
139. In the steady-state test case, the stress distribution is linear over the
depth (Equation 143); therefore, both the DSSi and the DSS§ give the exact solution
using one element over the vertical. The number of degrees of freedom (NDF) in the
120
finite element solution, (i.e., the number of simultaneous equations that must be
solved) is equal to the number of nodes used in the discretization (number of nodes =
1 + niimber of elements). The NDF in the spectral solution is equal to the number
of LPs used in the discretization. In both the DSS^ method and the DSS^ method,
the exact steady-€tate solution is obtained using two degrees of freedom.
140. Results from the periodic test case are shown in Figures 22 - 26. When
ft < 1, the DSS* is nearly exact using one finite element (two degrees of freedom)
over the depth, Figures 22a and 22b, 23a and 23b, 24a and 24b, and 25a and 25b.
For ft > 1, more than one finite element is required over the vertical for either DSS^
to converge to the analytical solution (Figures 22c and 22d, 23c and 23d, 24c and
24d, and 25c and 25d). Comparing these results to the DSS< results indicates that
both DSS^ methods require more degrees of freedom than the DSS» method to reach
the same level of convergence. The bottom stress plots presented in Figure 26
demonstrate the properties of the DSS* method further. In particular, they indicate
that the DSS^ is quite effective in the range ft ~ 10 or less. It may be possible to
extend this range to higher values of ft if an unequally spaced finite element grid is
used over the depth.
141. A time history of bottom stress for the transient test case is shown in
Figure 27. Comparing this to Figure 21a indicates that both DSS* methods require
four degrees of freedom to give a solution that is approximately equivalent to the
DSS‘ using three degrees of freedom.
142. In conclusion, new internal mode equations have been developed that
allow shear stress to be used as the dependent variable in the internal mode solution
and that yield a nearly tri-diagonal matrix system. While both DSS* require more
degrees of freedom than the DSS^ method to obtain comparable results for ft > 1,
(due to the use of linear finite elements in the DSS* versus spectral functions in the
DSSi), the matrix structure of the DSS* matrices makes this method much more
efficient than the DSS‘.
Implementation of Wave-Current Interaction in a DSS Model
143. It is often observed in lakes, coastal waters, and shelf waters that near
the bottom the orbital velocities associated with surface waves are as large as or
larger than the mean current velocity. In such cases the surface waves have a
significant effect on the bottom stress and the current profile. Several investigators
have developed theoretical models to account for this wave-current interaction. To
121
E*act Solution
Ld ^
u ^
si?
3
uj
3
Solution
124
I T,/ph 1 I r,/ph I
Figure 24. Vertical profiles of horizontal velocity magnitude for the periodic test case using the DSSi
Vertical profiles of horizontal velocity phase for the periodic test case using the DSS^
127
thoroughly assess the usefulness of the DSS approach, the effort required to implement
the Grant and Madsen (1979) model (GM model) with a DSS of the internal mode
equations has been considered. The GM model assumes that the mean current
velocity can be determined as follows:
a. inside the wave boundary layer, z < 6m,,
Ey = K jU^icwl 2 (200)
ju+cwl = ^ (201)
a no-slip boundary condition is applied at z = Zo, where Zq is the
physical bottom roughness
b. outside the wave boundary layer, z > S,,
Ey = « |U*c| z (202)
iu.c|=ijrcr= (203)
a no-slip boundary condition is applied at z = zoa, where Zoa is an
apparent bottom roughness experienced by the current due to the
wave-current interaction.
In these relations, ^ = 0.4 is the Von Kirman constant, Tc is the bottom stress due
to the current alone, Tw is the maximum wave-4nduced bottom stress during a wave
cycle, and is the thickness of the wave boundary layer.
144. The GM model can be included in a DSS of the internal mode equations
as follows.
g. Estimate Zoa and jU^d based on values at the previous time step.
b. Calculate Ey and use the DSS model to predict Tc-
c. Solve Equation 201 for lU*cw! using Tq from the previous step and
Tm, from Equation 53 in Grant and Madsen (1979). Since Tw is a
function of Equation 201 must be solved iteratively.
d. Determine Zoa using Equations 46 and 49 in Grant and Madsen
(1979).
g. Recalculate Ey using the new Tq. Use tliis and the new value of Zoa
in the DSS model to predict Tc- Go to step c. and iterate until Tq
(X)nverges.
145. Because two levels of iterations are required to implement the wave-
current interaction, it may be computationally infeasible to use this scheme in
128
practical model applications. It may be possible to simplify this procedure in two
ways in the proposed model. First, rather than iterate as described in step g., Zoa
can be calculated explicitly in time based only on results from the previous time step.
This should introduce little error into the solution if the time step is small enough
that changes in Zoa and Tq are relatively small. Second, following the suggestion of
Spaulding and Isaji (1987), Tw can be determined by neglecting the effect of the
current on the wave within the wave boundary layer. In this case
Tw/p = 0.5 fw lUbI Ub (204)
where fw is the wave friction factor (Jonsson and Carlsen 1976) and Ub is the
maximum bottom wave orbital velocity. For fully rough, turbulent flow, fw can be
determined from
— 1= + l08.» = ‘°6io ^ - 0.12 (205)
4 ^ fw 4 ^ fw
where Ab is the bottom excursion amplitude of the wave and kg is the Nikuradse
equivalent sand roughness of the bottom (typically Zq = ks/30).
146. The brief outline presented above suggests that the GM wave-current
interaction can easily be included in the DSS model. In fact, if the implementation
procedure outlined above for the DSS is compared with that described in Grant and
Madsen (1979) for a standard VS, it is evident that the DSS simplifies the use of the
GM model by eliminating the complications introduced by a quadratic slip bottom
boundary condition.
129
PART V: SUMMARY AND CONCLUSIONS
147. This report documents the theory and methodology behind the ADCIRC
(Advanced Circulation) model’s 2DDI (2-dimensional, depth-integrated) option and the
3DL (3-dimensional, local internal mode equation) option. ADCIRC is based on the
three-dimensional Reynold’s equations simplified using the hydrostatic pressure and the
Boussinesq approximations. Prior to their solution, the three-dimensional equations
are separated into a set of external mode equations (the two-dimensional, vertically
integrated equations) and a set of internal mode equations.
148. The external mode equations can be solved by themselves (the 2DDI
option) for depth-averaged velocity and free-surface elevation by parameterizing
bottom stress and momentum dispersion in terms of the depth-averaged velocity. Key
features of the external mode solution are the use of a generalized wave-continuity
equation (GWCE) formulation and the Galerkin-finite element (FE) method in space
using triangular or quadrilateral elements. The FE method provides maximum grid
flexibility and allows highly efficient numerical solutions to be obtained using model
domains that include complicated bathymetries and shoreline geometries that also
stretch considerable distances offshore to implement open-water boundary conditions.
Detailed analyses and testing of ADCIRC-2DDI have shown that it has good stability
characteristics, generates no spurious artificial modes, has minimal inherent numerical
damping, and efficiently separates the external mode equations into small systems of
algebraic equations with time-4ndependent matrices. Applications of the
ADCIRC-2DDI model to the English Channel and southern North Sea, the Gulf of
Mexico, Masonboro Inlet, and the New York Bight have shown that it is capable of
running month to year-long simulations while providing detailed intra-tidal
computations.
149. In stratified flows, Ekman layers, wind-driven flows in enclosed or semi-
enclosed basins, or flows affected by wave-current interaction in the boundary layer, it
is generally impossible to parameterize bottom stress and momentum dispersion in
terms of depth-averaged velocity. In such cases, it is necessary to solve the internal
mode equations for the vertical variation of horizontal velocity and use this to
evaluate the bottom stress and momentum dispersion terms in the external mode
equation. Due to the shallow water depths that characterize coastal and shelf
settings, the internal mode equations can often be simplified by dropping the
horizontal gradient terms. This gives internal mode equations that express the vertical
distribution of momentum at any horizontal position as a local balance between the
130
surface and bottom stresses, vertical momentum diffusion, the Coriolis ^orce, and local
inertia. The SDL model option is formulated using the simplified, local internal mode
equations. Existing numerical solutions of full or simplified internal mode equations
use velocity as the dependent variable. Therefore, it is necessary to use a fine
numerical discretization to resolve the sharp vertical gradients of velocity that occur
near the bottom boimdary and in wind-driven flows near the surface boundary.
During the course of the ADCIRC-3DL model development, a novel technique was
discovered that replaces velocity with shear stress as the dependent variable in the
internal mode equations. The resulting direct stress solution (DSS) allows physically
realistic boundary layers to be explicitly included in a three-dimensional model.
Detailed testing of the DSS method has demonstrated its considerable advantage over
standard velocity solutions and has led to an optimized DSS formulation. This
treatment of the internal mode equations should be invaluable for modeling coastal
and shelf circulation in which the bottom and surface boundary layers comprise a
significant portion of the water column and for modeling processes that are critically
dependent on boundary layer physics such as wave-current interaction, sediment
transport, oil spill movement, ice floe movement, energy dissipation, physical-biological
couplings, etc.
150. Considerable effort has gone into the development of ADCIRC to produce
a model that has simultaneous regional/local capabilities, as well as very high levels of
accuracy and efficiency. This has been achieved by combining extreme grid flexibility
with optimized formulations of the governing equations and numerical algorithms.
Together, these allow ADCIRC to run with improved physical realism and a
significant reduction in the computational cost of most presently existing circulation
models.
131
REFERENCES
Abbot, M. B. 1990. "Numerical Modeling for Coastal and Ocean Engineering,"
Handbook of Coastal and Ocean Engineering. J. Herbich, ed., Gulf Publishing.
ASCE. 1988a. "Turbulence Modeling of Surface Water Flow and Transport. Part I,"
J. Hydraulic Engineering. Vol 114, No. 9, pp 970-991.
_ . 1988b. "Turbulence Modeling of Surface Water Flow and Transport; Part II,"
J. Hydraulic Engineering. Vol 114, No. 9, pp 991-1014.
Baptista, A. M., Westerink, J. J., and Turner, P. J. 1989. "Tides in the English
Channel and Southern North Sea. A Frequency Domain Analysis Using Model
TEA-NL," Advances in Water Resources. Vol 12, No. 4, pp 166-183.
Bedford, K. W. 1984. "Selection of Turbulence and Mixing Parameterizations for
Estuary Models," Report to the US Army Engineer Waterways Experiment Station,
Environmental Impact Research Program, Vicksburg, MS.
Blumberg, A. F., and Mellor, G. L. 1987. "A Description of a Three-Dimensional
Coastal Ocean Circulation Model," Three-Dimensional Coastal Ocean Models. N.S.
Heaps, ed., AGU Press, Washington, DC, pp 1-16.
Bowden, K. F., Fairmairn, L. A., and Hughes, P. 1959. "The Distribution of Shearing
Stresses in a Tidal Current," Geophysical Journal of the Roval Astronomical Society.
Vol 2, pp 288-305.
Chu, W-S., Liou, J. -Y., and Flenniken, K. D. 1989. "Numerical Modebng of Tide
and Current in Central Puget Sound: Comparison of a Three-Dimensional and a
Depth-Averaged Model," Water Resources Research. Vol 25, No. 4, pp 721-734.
Connor, J. J, and Brebbia, C. A. 1976. Finite Element Techniques for Fluid Flow.
Newnes-Butterworths, London.
Csanady, G. T. 1976. "Mean Circulation in Shallow Seas," J. Geophysical Research.
Vol 81, No. 30, pp 5389-5399.
Davies, A. M. 1985. "Application of a Sigma Coordinate Sea Model to the Calculation
of Wind-Induced Currents," Continental Shelf Research. Vol 4, pp 389-423.
_ . 1987. "Spectral Models in Continental Shelf Oceanography,"
Three-Dimensional Coastal Ocean Models. N.S. Heaps, ed., AGU Press, Washington,
DC, pp 71-106.
_ . 1988. "On Formulating Two-Dimensional Vertically Integrated
Hydrodynamic Numerical Models with Enhanced Representation of Bed Stress,"
J. Geophysical Research. Vol 93, No. C2. pp 1241-1263.
_ . 1989. Modeling Marine Systems. Vol 1, CRC Press, Boca Raton, FL.
_ . 1991. "On the Accuracy of Finite Difference and Modal Methods for
Computing Tidal and Wind Wave Current Profiles," International J. Numerical
Methods in Fluids. Vol 12, No. 2, pp 101-124.
132
Davies, A. M., and Furnes, G. 1980. "Observed and Computed M2 Tidal Currents in
the North Sea,” J. Physical Oceanography. Vol 10, pp 237-257.
Davies, A. M., and Jones, J. E. 1990. "Application of a Three-Dimensional Turbulence
Energy Model to the Determination of Tidal Currents on the Northwest European
Shelf," J. Geophysical Research. Vol 95, No. CIO, pp 18,143-18,162.
Davies, A. M., and Jones, J. E. 1991. "On the Numerical Solution of the Turbulence
Energy Equations for Wave and Tidal Flows," International J. Numerical Methods in
Fluids. Vol 12, No. 1, pp 12^2.
Davies, A. M., and Owen, A. 1979. "Three-Dimensional Numerical Sea Model Ujing
the Galerkin Method with a Polynomial Basis Set," Applied Mathematical Modelling,
Vol 3, pp 421^28.
Davies, A. M., and Stephens, C. V. 1983. "Comparison of the Finite Difference and
Galerkin Methods as Applied to the Solution of the Hydrodynamic Equations,"
Applied Mathematical Modelling. Vol 7, pp 226-240.
Ferziger, J. H. 1987. "Simulation of Incompressible Turbulent Flows," J. Computa¬
tional Physics. Vol 69, pp 1-48.
Fischer, H. B., List, E. J., Koh, R. C. Y., Imberger, J., and Brooks, N. H. 1979.
"Mixing in Inland and Coastal Waters". Academic Press, New York.
Foreman, M. G. G. 1977. Manual for Tidal Heights Analysis and Prediction. Pacific
Marine Science Report 77-10, Institute of Ocean Sciences, Patricia Bay, Sidney, B.C
_ . 1983. "An Analysis of the Wave Equation Model for Finite
Element Tidal Computations," J. Computational Physics. Vol 52, pp 290-312.
Furnes, G. 1983. "A Three-Dimensional Numerical Sea Model with Eddy Viscosity
Varying Piecewise Linearly in the Vertical." Continental Shelf Research. Vol 2, No. 4,
pp 231-241.
Grant, W. D., and Madsen, O. S. 1979. "Combined Wave and Current Interaction
with a Rough Bottom," J. Geophysical Research. Vol 84, No. C4, pp 1797-1808.
_ . 1986. "The Continental Shelf Bottom Boundary
Layer," Annual Review of Fluid Mechanics. Vol 18, pp 265-305.
Gray, W. G. 1982. "Some Inadequacies of Finite Element Models as Simulators of
Two-Dimensional Circulation," Advances in Water Resources. Vol 5, pp 171-177.
_ . 1984. "On Normal Flow Boundary Conditions in Finite Element Codes
for Two-Dimensional Shallow Water Flow," International J. Numerical Methods in
Fluids. Vol 4, pp 99-104.
Gray, W. G., and Lynch, D. R. 1979. "On the Control of Noise in Finite Element
Tidal Computations: a Semi-Implidt Approach," Computers in Fluids. Vol 7,
pp 47-67.
133
Gresho, P. M., Lee, R. L., and Sani, R. L. 1987. "The Consistent Method for
Computing Derived Boundary Quantities When the Galerkin FEM is Used to Solve
Thermal and/or Fluids Problems," International J. Numerical Methods in Fluids.
Vol 7, pp 371-394.
Hansen, W. 1956. "Theorie zur Errechnung des Wasserstandes un der Strom ungen in
Randmeeren Nebst Anwendungen," Tellus. Vol 8.
Heaps, N. S. 1987. Three-Dimensional Coastal Ocean Models. AGU Press, Washington,
DC.
Hendershott, M. C. 1981. "Long Waves and Ocean Tides," Evolution of Physical
Oceanography. MIT Press, Cambridge, MA, pp 292-341.
Jenter, H. L., and Madsen, 0. S. 1989. "Bottom Stress in Wind-Driven
Depth-Averaged Coastal Flows," J. Physical Oceanography. Vol 19, No. 7,
pp 962-974.
Johns, B., and Oguz, T. 1987. "Turbulent Energy Closure Schemes," Three-
Dimensional Coastal Ocean Models. N.S. Heaps, ed., AGU Press, Washington, DC,
pp 17-40.
Jonsson, I. G., and Carlsen, N. A. 1976. "Experimental and Theoretical Investigations
in an Oscillatory Turbulent Boundary Layer," J. Hydraulic Research. Vol 14, No. 1,
pp 45-60.
Kinmark, I. 1985. "The Shallow Water Wave Equations: Formulation, Analysis and
Application," Lecture Notes in Engineering. Vol 15, C.A. Brebbia and S.A. Orszag,
eds., Springer-Verlag, Berlin.
Kolar, R. L., and Gray, W. G. 1990. "Shallow Water Modeling in Small Water
Bodies," Proceedings of the Eighth International Conference on Computational Methods
in Water Resources, pp 39-44.
Koutitas, C. 1987. "Three-Dimensional Models of Coastal Circulation: An Engineering
Viewpoint," Three-Dimensional Coastal Ocean Models. N.S. Heaps, ed., AGU Press,
WasMngton, DC, pp 107-124.
Launder, B. E. 1984. "Second-Moment Closure: Methodology and Practice," Simulation
of Turbulence Models and Their Applications. Vol 2, Collection de la Direction des
Etudes et Recherches, Electricite de France, editions Eyrolles, Paris.
Leendertse, J. J. 1967. "Aspects of a Computational Model for Long-Period Water
Wave Propagation," Memorandum RM5294-PR. Rand Corporation, Santa Monica, CA.
Luettich, R. A., Jr., and Westerink, J. J. 1991. "A Solution for the Vertical Variation
of Stress, ^ther Than Velocity, in a Three-Dimensional Circulation Model,"
International J. of Numerical Methods in Fluids, Vol 12, pp 911-928.
Lynch, D. R. 1978. "Finite Element Solution of the Shallow Water Equations," Ph.D.
dissertation, Department of Civil Engineering, Princeton University, Princeton, N.J.
134
Lynch, D. R. 1981. "Comparison of Spectral and Time-Stepping Approaches for Finite
Element Modeling of Tidd Circulation," Oceans 81 Conference Record. IEEE Publ.
No. 81CH1685-7, Boston, MA.
Lynch, D. R., and Gray, W. G. 1978. "Analytical Solutions for Computer Flow Model
Testing," J. of the Hydraulics Division. ASCE, No. HYIO, pp 1409-1428.
_ ^ _ . 1979. "A Wave Equation Model for Finite Element
Tidal Computations." Computers and Fluids. Vol 7, pp 207-228.
Lynch, D. R., and Officer, C. B. 1985. "Analytic Test Cases for Three-Dimensional
Hydrodynamic Models," International J. Numerical Methods in Fluids. Vol 5,
pp 529-^43.
Lynch, D. R., and Werner, F. E. 1987. "Three-Dimensional Hydrodynamics on Finite
Elements. Part I: Linearized Harmonic Model," International J. Numerical Methods in
Fluids. Vol 7, pp 871-909.
_ . 1991. "Three-Dimensional Hydrodynamics on Finite
Elements. Part II: Nonlinear Time-Stepping Model," International J. Numerical
Methods in Fluids. Vol 12, No. 6, pp 507-534.
Mellor, G. L., and Yamada, T. 1982. "Development of a Turbulence Closure Model for
Geophysical Fluid Problems," Review of Geophysics and Space Physics. Vol 20,
pp 851-875.
Monin, A., and Yaglom, A. 1971. Statistical Hydromechanics. Vol 1, MIT Press,
Cambridge, MA.
Nihoul, J. C. J., and Jamart, B. M. 1987. Three-Dimensional Models of Marine and
Estuarine Dynamics. Elsevier Science Publishing Co., Amsterdam.
Nihoul, J. C. J., and Djenidi, S. 1987. "Perspective in Three-Dimensional Modeling of
the Marine System," Three-Dimensional Models of Marine and Estuarine Dynamics.
J. J. Nihoul and B. M. Jamart, eds., Elsevier Science Publishing Co., Amsterdam,
pp 1-19.
Oey, L. -Y., Mellor, G. L., and Hires, R. I. 1985. "A Three-Dimensional Simulation
of the Hudson-Raritan Estuary. Part 1: Description of the Model and Modd
Simulations," J. Physical Oceanography. Vol 15, No. 12, pp 1676-1692.
Pinder, G. F., and Gray, W. G. 1977. Finite Element Simulation in Surface and
Subsurface Hydrology. Academic Press, New York.
Platzman, G. W. 1981. "Some Response Characteristics of Finite Element Tidal
Models," J. of Computational Physics. Vol 40, pp 36-63.
Poon, Y. -K. 1988. "A Two-Layer Coupled Hydrodynamic and Ice Floe Movement
Model", Sc.D. thesis, Department of Civil Engineering, Massachusetts Institute of
Technology, Cambridge, MA.
Reid, R. 0. 1990. "Water Level Changes - Tides and Storm Surges," Handbook of
Coastal and Ocean Engineering. Gulf Publishing Co.
135
Rodi, W. 1984. Turbulence Models and Thdr Application in Hydraulics - a State of
the Art Review. International Association for Hydraulic Research, Delft, Netherlands.
_ 1987. "Examples of Calculation Methods for Flow and Mixing in Stratibed
Fluids,*' J. Geophysical Research. Vol 92, No. C5, pp 5305-5328.
Schlichting, H. 1979. Boundary Laver Theory. 7th ed., McGraw-Hill, New York,
Schureman, P. 1941. Manual of Harmonic Analysis and Prediction of Tides. Special
Pub. No. 98, Coast and Geodetic Survey, US Dept, of Commerce, US Government
Printing Office, Washington, DC.
Schwiderski, E. W. 1980. "On Charting Global Ocean Tides," Review of Geophysics
and Space Physics. Vol 18, No. 1, pp 243-268.
Sheng, Y. P. 1983. "Mathematical Modeling of Three-Dimensional Coastal Currents
and Sediment Dispersion: Model Development and Application," Technical Report
CERC-83-2, US Army Engineer Waterways Experiment Station, Vicksburg, MS.
Sheng, Y. P., and Lick, W. 1980. "A Two-Mode Free-Surface Numerical Model for
the Three-Dimensional Time-Dependent Currents in Large Lakes," EPA Project
Report EPA-600/3-80-047.
Signell, R. P. 1989. "Tidal Dynamics and Dispersion Around Coastal Headlands,"
Report No. WHOI-89-38, WHOI-MIT Joint Program in Oceanography and
Oceanographic Engineering.
Simmons, T. J. 1974. "Verification of Numerical Models of Lake Ontario: Part I.
Circulation in Spring and Early Summer," J. Physical Oceanography. Vol 4,
pp 507^23.
Smith, L. H., and Cheng, R. T. 1987. "Tidal and Tidally Averaged Circulation
Characteristics of Suisan Bay, California," Water Resources Research. Vol 23, No. 1,
pp 143-155.
Spaulding, M. L. 1984. "A Vertically Averaged Circulation Model Using Boundary
Fitted Coordinates," J. Physical Oceanography. Vol 14, pp 973-982.
Spaulding, M. L., and Isaji, T. 1987. "Three Dimensional Continental Shelf
Hydrodynamics Model Including Wave Current Interaction," Three-Dimensional Models
of Marine and Estuarine Dynamics. J.J. Nihoul and B.M. Jamart, eds., Elsevier
Science Publishing Co., Amsterdam, pp 405-426.
Tee, K. -T. 1979. "The Structure of Three-Dimensional Tide-Generating Currents.
Part I: Oscillating Currents," J. Physical Oceanography. Vol 9, pp 930-^44.
_ 1987. "Simple Models to Simulate Three-Dimensional Tidal and Residual
Currents " Three-Dimensional Coastal Ocean Models. N.S. Heaps, ed., AGU Press,
Washington, DC, pp 125-148.
Vincent, P., and Le Provost, C. 1988. "Semidiurnal Tides in the Northeast Atlantic
from a Finite Element Numerical Model," J. Geophysical Research. Vol 93, No. Cl,
pp 543-555.
136
Walters, R. A. 1983. "Numerically Induced Oscillations in Finite Element
Approximations to the Shallow Water Equations," International J. of Numerical
Methods in Fluids. Vol 3, pp 591-604.
_ . 1984. "Finite Element Solution Methods for Circulation in Estuaries,"
Proceedings of the Fifth International Conference on Computational Methods in Water
Resources. J.P. Liable, et al., eds., Burlington, VT.
_ . 1987. "A Model for Tides and Currents in the English Channel and
Southern North Sea," Advances in Water Resources. Vol 10, pp 138-148.
Walters, R. A., and Carey, G. F. 1983. "Analysis of Spurious Oscillation Modes for
the Shallow Water and Navier-Stokes Equations," Computers and Fluids, Vol 11,
pp 51-68.
Wang, D. -P. 1982. "Development of a Three-Dimensional, Limited-Area (Island)
Shelf Circulation Model," J. Physical Oceanography. Vol 12, pp 605-617.
Wang, J. D., and Connor, J. J. 1975. "Mathematical Modeling of Near Coastal
Circulation," R.M. Parsons Laboratory Tech. Rep. 200, Massachusetts Institute of
Technology, Cambridge, MA.
Werner, F. E., and Lynch, D. R. 1987. "Field Verification of Wave Equation Tidal
Dynamics in the English Channel and Southern North Sea," Advances in Water
Resources. Vol 10, pp 184-193.
_ . 1988 (Feb). "Tides in the Southern North
Sea/English Channel; Data Files and Procedure for Reference Computations,"
Dartmouth College Numerical Laboratory Report, Dartmouth, NH.
_ _ ^ _ . 1989. "Harmonic Structure of English
Channel/Sourhern Bight Tides from a Wave Equation Simulation," Advances in Water
Resources. Vol 12, pp 121-142.
Westerink, J. J., Connor, J. J., and Stolzenbach, K. D. 1987. "A Primitive Pseudo
Wave Equation Formulation for Solving the Harmonic Shallow Water Equations,"
Advances in Water Resources. Vol 10, pp 188-199.
_ . 1988. "A Frequency-Time
Domain Finite Element Model for Tidal Circulation Based on the Least-Squares
Harmonic Analysis Method," International J. of Numerical Methods in Fluids. Vol 8,
pp 813-843.
Westerink, J. J., Stolzenbach, K. D., and Connor, J. J. 1989. "General Spectral
Computations of the Nonlinear Shallow Water Tidal Interactions within the Bight of
Abaco," J. Physical Oceanography. Vol 19, No. 9, pp 1348-1371.
Westerink, J. J., and Gray, W. G. 1991. "Progress in Surface Water Modeling,"
Reviews of Geophysics. Supplement, pp 210-219.
Westerink, J. J., Luettich, R. A., Jr., Baptista, A. M., Scheffner, N. W., and
Farrar, P. "Tide and Storm Surge Predictions Using a Finite Element Model,"
J. Hydraulic Engineering, in press.
137
Waterways Experiment Station Cataioging-in-Publication Data
Luettich, Richard A.
ADCIRC : an advanced three-dimensional circulation model for
shelves, coasts, and estuaries. Theory and methodology of ADCIRC-
2DDI and ADCIRC-3DL/ by R.A. Luettich, Jr., J.J. Westerink, Norman
W. Scheffner ; prepared for Department of the Army, US Army Corps of
Engineers.
141 p. : ill. ; 28 cm. — (Technical report ; DRP-92-6 rept. 1)
Includes bibliographical references.
1 . Ocean circulation — Mathematical models. 2. Hydrodynamics —
Mathematical models. 3. Ocean currents — Mathematical models. I.
Westerink, J. J. II. Scheffner, Norman W. III. United States. Army.
Corp® of Engineers. IV. U.S. Army Engineer Waterways Experiment
Station. V. Dredging Research Program. VI. Title. Vil. Title: ADCIRC:
an advanced three-dimensional circulation model for shelves, coasts,
and estuaries. Vlil. Trtle: Theory and methodology of ADCIRC-2DDi
and ADCIRC-3DL. IX. Series: Technical mport (U.S. Army Engineer
Waterways Experiment Station) ; DRP-92-6 rept.1 .
TA7 W34 no.DRP-92-6 rept.1
Live Music Archive
Librivox Free Audio
Metropolitan Museum
Cleveland Museum of Art
Internet Arcade
Console Living Room
Open Library
American Libraries
TV News
Understanding 9/11